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ABSTRACT 


Ongoing  interest  in  establishing  a  base  on  Mars  has  spurred  a  need  for  regular  and 
repeated  visits  to  the  red  planet  using  a  cycling  shuttle  to  transport  supplies,  equipment 
and  to  retrieve  surface  samples.  This  thesis  presents  an  approach  to  designing  an  optimal 
heliocentric  cycling  orbit,  or  cycler,  using  solar  sails.  Results  show  that  solar  sails  can  be 
used  to  significantly  reduce  Vj,  at  Mars  and  Earth.  For  example,  using  a  reasonably 
high  performance  solar  sail,  a  Vx  \Mar  =  2.5  km/s  is  possible  at  every  synodic  period  using 

a  two-dimensional  orbit  model.  Lower  performance  sails  were  also  modeled  resulting  in 
paths  that  behaved  more  like  a  ballistic  Aldrin  cycler  with  higher  V^s.  Double 
rendezvous  missions  were  explored  where  the  spacecraft  must  match  the  velocities  of 
both  Earth  and  Mars,  offering  promising  trajectories  for  Mars  sample  return  missions. 
The  solutions  to  these  missions,  although  not  necessarily  cyclers,  show  that  using  a  sail  to 
rendezvous  with  and  remain  near  Mars  for  an  optimal  amount  of  time  will  minimize  the 
total  transit  time  between  Earth  and  Mars.  General-purpose  dynamic  optimization 
software,  DIDO,  is  used  to  solve  the  optimal  control  problem  using  a  pseudospectral 
method  using  both  two-  and  three -dimensional  elliptic  orbit  models. 
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I.  INTRODUCTION 


Spacecraft  orbiting  along  cycling  trajectories  between  Earth  and  a  nearby  planet 
such  as  Mars  could  be  tremendously  valuable  for  missions  ranging  from  sample  returns 
to  ferrying  supplies  and  personnel  between  the  planets  like  escalators.  Naturally 
occurring  cycler  trajectories  would  be  ideal,  requiring  no  added  energy  other  than  that 
provided  by  the  gravity  fields  of  the  target  planets.  If  such  ‘free  ride”  paths  exist,  they 
could  host  spacecraft  making  an  infinite  number  of  round  trips  without  requiring 
propulsion  systems.  This  hypothetical  path  would  be  characterized  by  a  purely  natural 
Keplerian  trajectory  that  repeats  itself  endlessly  interrupted  only  by  planetary  swingbys. 
Such  cyclers  could  be  analyzed  for  various  target  planets,  however  Earth  -Mars  cyclers 
will  be  studied  in  this  thesis  because  of  relevance  to  establishing  a  base  on  our 
neighboring  red  planet.  To  this  end,  two  prominent  Earth-Mars  cycler  concepts  have 
been  proposed-  the  Aldrin  cycler1  and  VISIT  cycler2.  The  Aldrin  cycler  uses  gravity 
assists  from  both  planets  and  small  well-timed  A  Vs  to  maintain  a  continuous  cycling 
trajectory.  While  the  revisit  times  are  appealing  (7  round  trips  in  15  years),  on-board 
propellant  is  required  to  provide  a  15  year  cumulative  AV  of  approximately  2  km/s. 
Because  no  fuel  is  expended  trying  to  reduce  speed  near  a  planet,  hyperbolic  excess 
speeds  are  high,  in  excess  of  5  km/s.  A  VISIT  cycler,  on  the  other  hand,  uses  neither 
gravity  assists  nor  fuel  for  AV  bums  for  orbit  shaping,  but  rather  resides  in  a  natural 
heliocentric  orbit  that  makes  regular  passes  of  both  Earth  and  Mars.  The  advantage  to 
this  cycler  is  that  no  propellant  is  required  to  maintain  the  orbit  for  up  to  20  years,  but  the 
primary  disadvantages  are  the  rather  long  revisit  times  (3  Earth  visits  and  4  Mars  visits  in 
15  years)  and  large  approach  distances  needed  to  avoid  unwanted  gravitational 
interaction.  The  purpose  of  this  paper  is  to  present  an  Earth -Mars  cycler  concept  using 
solar  sails  that  has  the  advantages  of  both  cyclers  without  the  on  -board  propellant. 

Because  it  has  been  shown  that  a  nearly  ballistic  cycler  orbit  can  be  maintained 
for  15  years',  it  is  reasonable  to  speculate  that  a  solar  sail  cycler  is  not  only  feasible,  but 
more  capable  since  it  provides  free  controls.  With  solar  radiation  pressure  provid  ing  free 
thrust  to  the  solar  sail,  a  cycler  orbit  may  be  maintained  with  short  revisit  times  as  well  as 


1 


slow  planetary  approach  speeds.  Small  approach  velocities  to  Earth  and  Mars  at  short 
distances  are  desirable  for  operational  reasons  and  because  they  exhibit  attractive  orbit - 
shaping  characteristics  due  to  large  angle  swingby  deflections. 

Given  that  a  solar  sail  requires  no  propellant,  we  are  free  to  determine  a  cost 
function  that  does  not  minimize  the  fuel.  What  then  defines  an  optimal  solar  sail  cycle  r? 
Often  the  time  of  flight  or  an  undesirable  cycler  characteristic  like  large  approach 
velocities  may  be  minimized.  A  design  parameter  of  the  sail  itself,  such  as  the  size  of  the 
sail  may  be  the  set  cost  function.  This  thesis  presents  a  framework  for  solving  such 
optimal  low-thrust  cycler  trajectories  using  the  general-purpose  dynamic  optimization 
software  DIDO3.  It  is  worth  mentioning  that  although  solar  sails  are  chosen  as  the 
propulsion  method  and  cycling  orbits  are  selected  as  the  mission  for  this  thesis,  the 
framework  presented  here  is  applicable  to  any  type  of  propulsion  system  and  mission. 
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II.  APPROACH  TO  OPTIMAL  CONTROL  SOLUTIONS 


A.  ROADMAP 

At  the  outset  of  this  thesis  work,  there  were  no  known  works  related  to  optimal 
Earth-Mars  solar  sail  cyclers,  therefore  a  cautious  approach  to  the  optimal  control 
solution  was  necessary.  Starting  with  “simple”  two  point  boundary  value  astrodynamic 
problems,  a  systematic  approach  was  used  that  gradually  increased  the  model  complexity 
to  arrive  at  solutions.  These  solutions  were  used  to  both  validate  the  numerical 
optimization  method  and  to  obtain  milestone  profiles  providing  initial  guesses  and 
bounds  to  more  complex  problems.  Benchmark  optimal  control  problems  were  solved 
using  DIDO,  and  the  output  solutions  were  compared  with  known  solutions .  Benchmark 
problems  included  Mars  flyby.  Mars  rendezvous  and  Earth-Mars-Earth  double 
rendezvous  missions.  Numerical  results  were  also  propagated  using  initial  conditions  and 
the  output  control  histories  for  further  validation.  For  many  of  the  simpler  problems  (i.e. 
no  intermediate  event  conditions)  a  costate  profile  was  obtained  corresponding  to  all  the 
Lagrange  multipliers  at  certain  points  in  time  along  the  path.  Using  the  costates  with  the 
derived  solar  sail  steering  control  law  enabled  a  check  of  compatibility  between  the 
DIDO  costate  and  output  control  histories.  The  resulting  Hamiltonian  output  was  also 
available  to  check  for  optimality.  After  the  benchmark  problem  solutions  had  been 
verified,  additional  constraints  and  event  conditions  were  inc  luded  to  obtain  the  more 
complex  cycler  trajectory  that  includes  planetary  swingbys  of  Earth  and  Mars. 

The  above  approach  was  used  for  both  two-dimensional  and  three- dimensional 
models.  The  two-dimensional  solutions  approximated  both  Earth  and  Mars  orbits  as 
being  circular  and  coplanar  serving  as  an  initial  estimate  for  the  higher  fidelity  models. 
Fortunately,  the  relative  inclination  of  Mars’  orbit  with  respect  to  the  Earth  ecliptic  is 
only  1.85°  and  both  orbit  eccentricities  are  not  too  large,  so  this  approximation  is  not  too 
bad.  Increasing  the  fidelity  of  the  cycler  model,  Mars  was  given  an  elliptic,  but  coplanar 
orbit  to  observe  the  differences  in  the  optimal  paths.  Finally,  Earth  and  Mars  orbits  were 
inclined  and  both  eccentricities  were  taken  into  account.  The  results  from  the  battery  of 
2-D  problems  provided  bounds  and  initial  guesses  for  their  3-D  counterpart  problems. 
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Comparison  of  the  2-D  and  3-D  solutions  to  the  same  optimal  control  problem  provided 
some  interesting  insights  into  astrodynamic  optimization. 

1.  Coordinate  Systems 

To  improve  the  performance  of  the  numerical  analysis,  it  served  well  to  choose  a 
coordinate  system  with  states  that  change  slowly  with  time  and  dynamics  that  have  few 
or  no  singularities.  Using  a  coordinate  system  in  which  states  vary  relatively  slowly  with 
time  seemed  to  allow  accurate  solutions  with  fewer  optimization  “node”  points,  or  points 
along  the  path  that  are  used  to  approximate  the  continuous  states  and  controls  in  the 
numeric  optimization  process.  Additionally,  singularity  avoidance  is  necessary  to  obtain 
feasible  solutions. 


For  the  two-dimensional  model,  a  heliocentric  polar  coordinate  system  was  used 
as  shown  in  Figure  1 .  Because  the  maximum  solar  sail  thrust  is  inversely  proportional  to 
the  distance  from  the  sun,  this  coordinate  system  was  well  suited  for  this  problem.  The 
only  singularity  in  the  polar  equations  of  motion  occurs  at  the  center  of  the  sun  -  an  easily 
avoidable  situation  with  proper  state  bounds  (see  Appendix  A).  Velocity  components 
were  expressed  in  terms  of  radial  velocity  and  along -track  velocity  (ref  4  p.  43).  The 
along- track  velocity  will  be  referred  to  as  the  transverse  velocity  and  is  defined  as  being 
normal  to  radial  position  vector  in  the  orbital  plane.  Since  the  target  orbits  were 
approximated  as  circular  and  coplanar  in  the  2-D  model,  there  was  no  need  to  define  a 
principle  direction  (i.e.  along  vernal  equinox). 
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Likewise,  for  the  3-D  model  spherical  coordinates  were  chosen  to  represent  3-D 
states  and  dynamics  because  the  resulting  optimal  state  and  control  profiles  were 
expected  to  vary  slowly  with  time  allowing  excellent  numerical  approximations  using 
relatively  few  nodes.  Furthermore,  these  coordinates  made  it  very  simple  to  model  the 
thrust  due  to  the  sail  since  thrust  is  dependent  on  radial  distance,  r  (see  Appendix  B).  By 
comparison,  states  of  orbiting  bodies  in  Cartesian  coordinates  oscillate.  This  could  have 
lead  to  an  undesirable  condition  in  case  there  turned  out  to  be  too  many  cycles  in  the 
solution  for  the  number  of  nodes  to  approximate  accurately.  It  should  be  kept  in  mind, 
however,  that  Cartesian  coordinates  present  perhaps  the  best  choice  when  more 
perturbations  are  desired  in  the  dynamics  model  (in  this  case,  more  nodes  should  be 
used).  Equinoctial  elements  provide  a  set  of  singularity -free  slowly  changing  variables, 
but  because  it  was  desirable  to  keep  boundary  conditions  simple  and  “intuitive”,  they 
were  not  used  (ref  4  p.  142). 

2.  Scaling 

With  the  coordinate  system  chosen  to  reduce  the  number  of  optimization  nodes 
required  to  accurately  approximate  states  and  controls,  it  became  imperative  to  select 
units  that  would  reduce  numerical  conditioning  by  avoiding  the  extremely  large  and  small 
numbers.  There  is  no  need  to  confine  the  states  to  familiar  units  such  as  miles, 
kilometers,  seconds,  etc.  A  useful  unit  system  for  modeling  interplanetary  travel  is  the 
canonical  system.  Choosing  a  distance  unit  (DU)  to  be  1  A.U.  and  the  normalized  solar 
gravitational  parameter  (X  to  be  1,  other  state  units  may  be  calibrated  accordingly.  All 
states  and  controls  have  units  that  are  of  the  same  order  of  magnitude  and  dimensio  ns  are 
scaled  by  a  factor  of  some  characteristic  dimension.  This  non-dimensionalization  process 
also  enhances  visualizing  metric  relationships.  For  instance,  the  circular  Earth  orbit 
radius  is  1  DU  while  the  Mars  orbit  radius  is  1.524  DUs,  over  50%  greater  than  Earth’s. 
Heliocentric  canonical  units  are  summarized  in  Table  1. 
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Dimension 

Non-dimensional  unit 

inks  equivalent 

Comments 

Distance  unit 

DU 

1.496E8  km 

Semi-major  axis  of  Earth  orbit  =  1 
AU. 

Time  unit 

TU 

5.023E6  sec 

Time  to  traverse  1  radian  along 
circular  orbit  at  1  AU 

Speed 

VU  =  DU/TU 

29.78  km/s 

Velocity  of  body  in  1  AU  circular 
orbit. 

Gravitational 

parameter 

Solar  gravitation  parameter 

Table  1  Canonical  Units  Used  to  Non-Dimensionalize  the  States  and 

Parameters. 

3.  Bounding  the  Problems 

Before  approaching  the  full  cycler  model  accounting  for  gravity  assists,  simpler 
two-point  boundary  value  problems  were  analyzed.  As  already  mentioned,  this  permits 
validation  of  the  optimization  software.  Additionally,  important  insights  may  be  gained 
by  observing  the  behavior  of  simpler  benchmark  problems.  For  both  the  2-D  and  3-D 
models,  three  benchmark  problems  were  solved  prior  to  the  coding  of  the  cycler  models. 
These  are  the  Mars  flyby.  Mars  rendezvous,  and  Earth -Mars -Earth  double  rendezvous 
round  trip  missions.  The  Mars  flyby  mission  is  a  simple  minimum  time  trajectory  from 
Earth  orbit  to  Mars  orbit.  This  solution  provides  the  absolute  shortest  time  that  a  solar 
sail  of  a  given  performance  can  reach  the  Mars  orbit.  Other  missions  that  include  an 
Earth  to  Mars  transit  may  use  this  time  as  a  lower  bound  since  it  is  impossible  to  reach 
Mars  in  less  time.  A  Mars  rendezvous  mission  is  one  in  which  the  solar  sail  is  tasked  to 
reach  Mars  at  its  local  orbital  velocity  (V„  =  0  with  no  Mars  gravity  acting)  in  minimum 
time.  The  resulting  profile  is  not  only  useful  in  observing  optimal  rendezvous  behavior, 
but  also  for  bounding  the  double  rendezvous  problem  that  followed.  The  double 
rendezvous  mission  required  the  spacecraft  to  rendezvous  with  Mars,  then  track  the 
planet  for  a  time  (with  no  effect  due  to  Mars  gravity),  and  then  return  to  rendezvous  with 
Earth.  This  mission  introduced  phasing  constraints  to  ensure  that  Earth  was  there  when 
the  spacecraft  returned.  With  the  double  rendezvous  bounds  established  in  this  manner,  a 
state  discontinuity  was  introduced  representing  a  gravity  assist  from  the  Mars  encounter 
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making  a  zero-sphere  patched  conic  cycler  model  complete.  At  this  point,  all  the  models 
and  conditions  had  been  established  for  a  solar  sail  cycler  and  coded  for  use  with  the 
numeric  optimization  software,  DIDO. 

4.  Numerical  Analysis,  terminology  and  DIDO  3 

Ideally,  the  solution  to  the  optimal  trajectory  would  view  the  entire  space -time 
path  of  a  spacecraft  at  once  and  simultaneously  locate  the  optimal  position  for  each  point 
along  that  path  meeting  certain  boundary  conditions  and  dynamic  constraints  (Figure  2 
and  Figure  3).  Limiting  ourselves  to  current  computer  methods,  however,  we  settle  for 
using  certain  points  along  the  space  time  path  that  will  enable  us  to  approximate  the 
continuous  path,  thus  dscretizing  the  trajectory  into  a  workable  form  for  a  non-linear 
programming  (NLP)  solver.  An  Euler  method  would  provide  an  approximation  that 
would  be  burdened  by  large  errors.  It  was  desired  for  this  complex  problem  to  use  the 
best  numerical  approximations  possible  in  solving  the  optimal  control  problem  (OCP). 


Figure  2  Path  and  Target  manifolds  (orbit  positions) 
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Figure  3  Path  Meeting  Event  Manifolds  (Planet  Encounters)  in  Space  (x-y)  and 

Time. 

The  numerical  solution  to  the  optimal  control  problem  presented  in  this  thesis 
uses  the  software  DIDO  to  interface  with  a  generic  third  party  NLP  solver.  The  NLP 
solver  requires  an  NLP  variable  vector,  constraints  and  a  cost  function.  DIDO  provides 
these  requirements  set  up  in  a  numerically  optimal  fashion.  All  state  and  control  path 
information  are  captured  at  unevenly  distributed  node  points,  thus  the  entire  dynamic 
problem  can  be  viewed  as  a  static  problem  with  the  complete  trajectory  history  at  selected 
node  points  manifested  as  a  single  variable  vector.  The  node  points  are  selected  in  such  a 
manner  that  the  state  and  control  approximation  error  is  minimized.  These  numerically 
optimal  node  points  are  called  Legendre-Gauss-Labatto  (LGL)  points  and  correspond  to 
the  zeros  of  the  first  derivative  of  the  Nth  degree  Legendre  polynomial'’'6.  The 
continuous  states  and  controls  are  approximated  by  using  the  variables  at  the 
predetermined  LGL  node  points  and  Lagrange  interpolating  polynomials,  which  are  well 
suited  for  interpolation  between  unevenly  distributed  data  points. 
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Dynamic  constraints  are  imposed  by  forcing  the  derivatives  of  the  approximated 
states  to  equal  user-defined  differential  equations  (equations  of  motion)  at  all  the  node 
points.  Obtaining  the  derivative  of  the  state  variable  vector  is  accomplished  using  a 
differentiation  matrix  defined  in  ref  5.  This  process  of  discretization  and  differentiation 
is  called  the  Legendre  pseudospectral  method  and  has  been  used  successfully  in  a  variety 
of  applications  (refs  5-8). 

Initial,  intermediate  and  final  event  conditions  a'e  supplied  by  the  user  and 
structured  into  a  boundary  condition  vector  that  DIDO  provides  as  constraints  to  a  third  - 
party  NLP  solver.  The  approach  used  in  this  numerical  solution  to  optimal  control  does 
not  require  propagation  or  shooting  type  methods.  An  initial  guess  is  required,  however 
it  does  not  need  to  be  a  good  one,  or  even  a  feasible  one.  Typically  the  end  conditions 
are  supplied  and  DIDO  interpolates  between  these  points  as  a  guess. 

To  assist  in  the  reading  of  the  numerical  methods  used  in  this  thesis,  some 
definitions  from  ref  3  are  supplied  below. 

Definitions: 

Nodes  are  discrete  points  along  the  path  where  states  and  controls  are  defined. 
The  entire  continuous  path  is  represented  by  the  states  and  controls  at  th  ese  discrete 
points.  The  states  and  controls  at  the  all  the  node  points  form  part  of  a  vector  making  up 
the  NLP  variable  that  inputs  to  the  NLP  solver.  Node  locations  are  unevenly  distributed 
and  are  pre-determined  at  the  LGL  points.  LGL  points  occur  at  the  zeros  of  the  first 
derivative  of  the  Legendre  polynomials  in  order  to  minimize  the  mean -squared  error  in 
Lagrange  interpolating  polynomial  approximation  of  states  and  controls. 

Knots  are  double  node  points  that  are  part  of  the  optimization  process.  T  rajectory 
end  points  are  also  called  knots.  Interior  knots  are  used  to  define  intermediate  events 
where  there  may  or  may  not  be  state  discontinuities.  LGL  node  distribution  always 
concentrates  more  nodes  next  to  knots.  Intermediate  or  interior  knots  have  a  left  and 
right  side  where  state  constraints  are  imposed. 

Events  are  the  set  of  conditions  that  define  the  manifolds  which  constrain  the  path 

(at  the  knots).  The  event  manifolds  usually  consist  of  the  boundary  conditions;  in  this 
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thesis  they  represent  the  locations  and  velocities  of  Earth  and  Mars  along  their  respective 
orbits  in  space  and  time.  Within  the  DIDO  structure,  knots  are  the  primary  way  to 
represent  an  event  condition.  An  example  of  event  boundary  conditions  is  shown  below 
where  e  represents  the  event  vector,  and  ei  and  eu  represent  the  lower  and  upper  bounds  of 
the  vector  respectively. 


such  that 

e;  fCelx^Xxlf.kxO^p,^/,  (f  )<e„ 

where  the  vector  p  contains  any  static  parameters  that  may  constrain  the  events.  Events 
form  the  constraints  that  DIDO  supplies  to  the  generic  NLP  solver  and  are  listed  in  tables 
for  the  various  OCPs  in  this  thesis. 

Path  constraints ,  g ,  are  the  state -control  constraints  that  couple  control  and  state 
variables  in  the  following  manner. 

g,  <g(x(f),u;f),f)<gH 

For  numerical  reasons  it  is  sometimes  advantageous  to  represent  the  solar  sail 
controls  as  the  sine  and  cosine  components  of  the  cone  and  clock  control  angles  instead 
of  the  actual  angles  themselves.  The  path  function  in  this  case  may  be  represented  as 

ll .  -T  u 0  -  —  1 

g  =  where  g,  =  g  =  0  forming  an  equality  constraint.  The  actual  control 

U{  +  It.,'  —  1 

angles  may  be  recovered  by  taking  the  arctangent  of  coupled  controls  (using  Matlab’s 
atan2  function  serves  well  here). 

Bounds  are  the  upper  and  lower  limits  on  states,  controls  or  times.  The  resulting 
restriction  is  in  the  form  of  an  equality  or  inequality  constraint.  For  instance,  a  set  of 
state  variables  may  be  bounded  by  x  e  S  c  R ,v  usually  written  in  the  form  x  e  [x  lh  ,x  ub  \ 

or  xlb  <  x  <  xuh .  Equality  constraints  are  implemented  simply  by  making  xlh  =  xub . 
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Within  DIDO,  states  and  controls  are  discretized,  so  variables  are  bounded  at  the  node 
points. 


Event  Constraints  are  those  linear  or  non-linear  conditions  that  restrict  states  at 
the  events.  For  instance,  in  an  orbital  rendezvous  problem,  the  final  spacecraft  velocity  is 


coupled  with  the  final  radius  of  the  target  planet  according  to  the  relation 


forming  a  non-linear  event  constraint.  When  a  state  at  a  particular  event  condition  is 
unconstrained,  it  is  considered  free.  For  example,  when  intercepting  Mars  in  minimum 
time,  we  do  not  specify  where  in  space  or  time  the  intercept  must  occur,  thus  final  time 
and  final  angular  displacement  are  free  variables. 


B.  VALIDATING  SOLUTT  ONS 

Once  a  solution  to  a  particular  trajectory  was  obtained,  it  was  important  to  see  if 
the  path  was  physically  feasible  and  if  so  was  it  optimal.  Checking  for  feasibility  proved 
to  be  relatively  straightforward,  while  verifying  optimality  was  a  bit  thornier.  This 
process  linked  the  stark  numerical  optimization  solution  with  a  physical  spacecraft 
trajectory  that  “made  sense”. 

1.  Feasibility 

A  glance  at  the  state  data  at  the  knots  would  verify  that  the  event  constraints  had 
been  met.  In  order  to  verify  that  the  entire  state  history  conformed  to  physical  laws 
governed  by  the  equations  of  motion  we  use  a  propagator.  Taking  the  control  history 
generated  by  the  numeric  optimization  and  the  initial  conditions,  the  path  of  the  sail  could 
be  propagated  by  means  of  a  numeric  ordinary  differential  equation  solver  on  the  non¬ 
linear  equations  of  motion.  The  propagation  used  the  same  dynamics  employed  by  the 
DIDO  code  to  verify  that  DIDO  was  properly  applying  the  dynamic  constraints  to  the 
path.  To  verify  the  equations  of  motion  in  the  dynamic  model  us  ed  by  DIDO,  the 
equations  were  propagated  without  controls  to  generate  familiar  Keplarian  orbits.  The 
state  and  control  history  figures  in  this  thesis  generally  show  the  DIDO  states  represented 
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by  markers  at  the  discrete  node  points  and  the  propagated  solution  as  a  line  that  will  pass 
through  DIDO’s  markers  when  the  solutions  are  consistent. 


2.  Optimality 

After  testing  for  feasibility,  the  solution  was  checked  for  optimality.  The  primary 
method  of  verifying  that  an  optimal  trajectory  had  been  achieved  was  observing  the 
behavior  of  the  Hamiltonian,  H(x,  u,  X ,  t).  DIDO-derived  Hamiltonian  values  were 
generated  at  the  predetermined  node  points  for  simple  problems  (i.e.  no  intermediate 
knots)  and  were  used  to  verify  minimum  time  optimal  controls.  Since  none  of  the 
problems  posed  in  this  thesis  have  a  Hamiltonian  as  explicit  function  of  time,  then  we  can 
say 

dH__dH__  o 
dt  dt 


Hamiltonians  that  are  constant  with  respect  to  time  are  necessary  (but  not 
sufficient)  conditions  for  optimality  in  these  cases.  Moreover,  the  optimal  control 
solution  is  the  one  where  H  is  minimized  with  respect  to  controls,  u  thus  defining  an  NLP 
problem  with  the  Lagrangian  of  the  Hamiltonian  (for  details  see  Appendix  A).  For 
minimum  time  problems,  H  =  —1  for  all  time.  The  first  and  second  order  necessary  (but 
not  sufficient)  conditions  for  optimality  using  the  Hamiltonian  are  written  as 


,«I  ,')LI 


Additionally,  dual  variable  outputs  were  used  with  the  derived  solar  sail  control 
law  to  produce  controls  that  could  be  compared  with  the  DIDO-derived  controls.  This 
would  test  if  the  optimal  controls  were  conforming  to  the  analytical  optimal  control 
steering  law  solution.  Occasionally  there  existed  one  or  more  local  minimums  in  which 
case  varying  the  initial  guess  assisted  in  finding  a  better  local  minimum. 


3.  Comparison  with  Other  Optimal  Solutions 

Often  it  serves  well  to  verify  DIDO-generated  optimal  solutions  that  use  the 
pseudos  pectral  method  with  solutions  created  using  other  methods.  As  mentioned  before, 
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there  are  no  known  optimal  control  solutions  to  the  actual  solar  sail  cycler,  however  there 
are  solutions  to  the  simpler  “benchmark”  problems.  One  of  the  reasons  for  using  the 
methodical  build-up  approach  using  simple  building  blocks  in  the  form  of  benchmark 
problems  was  to  verify  that  the  coded  models  were  producing  verifiable  results  before 
proceeding  to  the  uncharted  waters.  Comparisons  are  presented  in  the  benchmark 
problems  when  other  solutions  are  available. 
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III.  DEVELOPING  THE  OPTIMAL  TWO-DIMENSIONAL 

CYCLER 


A.  TWO-DIMENSIONAL  MODELS 


1.  Sail  Model 

A  solar  sail  is  a  very  large  thin  lightweight  structure  that  transfers  momentum 
from  inbound  solar  photons  to  the  spacecraft.  Newton’s  laws  dictate  that  the  sail  will 
accelerate  twice  as  fast  when  photons  are  reflected  than  when  they  are  absorbed.  Thus, 
for  maximum  effectiveness  the  sail  must  be  made  of  a  very  reflective  material.  For  this 
analysis,  the  sail  is  modeled  as  a  flat  film  (non -billowing  sail)  with  a  perfectly  specular - 
reflecting  surface  on  one  side  only.  To  characterize  sail  performance,  we  use  the  sail 
lightness  number,  (3  ,  defined  as  the  ratio  of  the  acceleration  from  solar  radiation  pressure 
(SRP)  to  the  acceleration  from  the  sun’s  gravity  (ref  9  pp.  38-40). 


P=- 


2  WR 


me 


■Acosta 


where  We  is  the  solar  energy  flux  at  1  AU,  R  is  the  Sun-Earth  distance,  A  is  the  sail 
area,  oc  is  the  angle  of  the  sail  with  respect  to  the  Sun  and  p  is  the  solar  gravitational 
parameter9.  Both  accelerations  are  proportional  to  the  inverse  square  of  the  radial 
distance  from  the  sun,  so  (3  is  an  apt  indicator  of  sail  performance  independent  of 
location.  A  lightness  number  of  (3  =.17  was  used  to  model  a  reasonably  high 
performance  solar  sail  for  most  problems  in  this  thesis  for  comparison  to  other  solutions 
using  the  same  sail10.  Often  the  characteristic  acceleration11,  aQ,  is  used  to  describe  the 
performance  of  the  sail  instead  the  sail  lightness  number.  This  characteristic  acceleration 
is  the  change  in  velocity  that  the  sail  would  experience  at  a  one  AU  distance  from  the  sun 
when  the  sail  exposes  all  of  its  area  toward  the  sun.  The  relationship  between  the 
lightness  number  and  the  characteristic  acceleration  is 
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=fl0 


For  example,  a  sail  lightness  number  of (3  =.17  corresponds  to  a  characteristic 

YYlYYl 

acceleration  of  a0  =  1  — — . 


The  sail  control  angle,  oc ,  is  defined  as  the  angle  between  the  sun- sail  position 
vector,  r,  and  the  sail  normal,  n  called  the  cone  angle  (Figure  4).  This  angle  determines 
both  the  magnitude  and  direction  of  a  solar  force  imparted  to  the  spacecraft.  Notice  that 

when  a  =  0 ,  maximum  thrust  is  achieved  and  when  a  =  90° ,  thrust  is  zero.  The  sail  can 
never  impart  a  force  toward  the  sun  (only  gravity  can  do  this).  To  ensure  the  sail  will  not 
flip  over  exposing  its  non- reflective  side,  the  control  angle  is  bounded  by 


Figure  4  Sail  normal,  solar  radial  vector  and  control  angle  a. 

2.  Dynamic  Model 

Using  the  polar  coordinates  in  the  heliocentric  frame,  the  two-dimensional 
equations  of  motion  are  derived  in  Appendix  R  Given  the  state  vector  x  =  [r,0  ,u,v]T , 
the  equations  of  motion  are  expressed  as  x  =  f  (x  ,u )  where  control  u  =  a  , 


16 


(2) 


u 

v 


V 

r 


ix  Bjx  3 

■ cos  a 
r  r~ 


uv  Bix  , 

- +  — ^-cos  a  sina 

r  r~ 


The  states  r,  0,  u,  and  v  are  the  radial  distance,  angular  displacement,  radial 
velocity  and  transverse  velocity  respectively.  The  angular  displacement,  0,  is  decoupled 
from  the  other  equations  and  may  be  excluded  in  special  cases  making  x  e  R:l .  For 
numerical  analysis  of  some  problems,  0  was  retained  as  a  state  for  proper  intercept 
phasing  during  encounter  events  (more  on  that  later). 

Having  thus  imposed  the  dynamic  constraints  between  the  desired  manifolds,  we 
seek  to  constrain  the  conditions  at  the  manifolds.  This  is  accomplished  by  setting  the 
constraints  of  the  desired  events. 


3.  Events  Model 

All  the  benchmark  missions  modeled  here  share  the  same  cost  function,  dynamic 
equations  of  motion  and  sail  control  model.  The  distinguishing  feature  of  each  mission  is 
its  particular  event  function.  The  events  are  arranged  into  an  event  function  vector  e  sue  h 
that  ee  R Np  where  N p  is  the  number  of  event  constraints.  Elements  of  e  may  be  linear 

or  nonlinear  expressions  and  may  be  bounded  as  equality  or  inequality  constraints.  These 
event  constraints  for  the  flyby  and  rendezvous  missions  are  summarized  in  the  following 
table  where  am  is  the  semi-major  axis  of  Mars’  orbit  (1.524  AU). 
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Event 

Flyby 

[lb,ub] 

Rendezvous 

[lb,ub] 

ro 

[U] 

[1,1] 

% 

[0,0] 

[0,0] 

uo 

[0,0] 

[0,0] 

VQ 

[1,1] 

[1,1] 

rf 

6/ 

free 

free 

Uf 

free 

0 

'V~7T 

free 

[0,0] 

Table  2  Event  Conditions  for  Flyby  and  Rendezvous  Problems 

Notice  that  all  the  initial  states  (denoted  by  the  ‘O’  subscript)  are  equality 
constraints  with  the  same  lower  and  upper  bounds  ([ lb,ub ])  indicating  that  the  spacecraft 
will  launch  from  a  circular  Earth  orbit  at  Earth’s  orbital  speed.  The  final  rendezvous 
along- track  velocity  (denoted  with  the  “f’  subscript)  is  the  only  non-linear  constraint. 
These  events  are  assembled  within  DIDO  into  a  c  onstraint  vector, 
e  =  [r0  0O  m(i  ,v0  ff  Qf,uf,v/rf  - 1 17 .  The  last  element  is  written  in  a  different  form  to 
avoid  the  +/  -  ambiguity  using  the  square  root. 


B.  THE  OPTIMAL  CONTROL  PROBLEM 

We  now  come  to  posing  the  solar  sail  optimal  control  problem  (OCP).  This  is 
where  one  must  question  what  exactly  are  we  trying  to  do.  With  most  conventional 
propulsion  schemes,  it  is  desirable  to  minimize  fuel  consumption.  Solar  sails,  however, 
require  no  propellant  so  the  goal  of  most  of  the  solar  sail  trajectories  presented  her  e  is  to 
accomplish  a  given  mission  in  minimum  time.  The  general  OCP  is  stated  as  follows3. 

Minimize  the  Bolza  cost  functional: 
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/[xQ, uQ,  t0,t  ■  pi  =  E(x(t0),x(tf ), t0,  t  ■  p)  +  f  E(x(t), u (f), t  ;p  )dt 


Subject  to: 

Dynamic  constraints  x  =  f  (x(f),u(/) ) 

where  f  is  given  in  equation  (2)  and  u  are  the  controls, 

Path  constraints  g,  <  g(x(f),u(  f),  t)  <  gu 

Event  constraints  e,  <e(x0,  x.,xp  ta  f ,  tj,  p)  <eu 

and  bounds  on  state  and  control  variables 

x;<x(r)<x„ 
u;  < u(f) <U( 

With  an  established  coordinate  system  and  scaling  system,  we  are  postured  to 
model  the  sail,  dynamics  and  events.  These  models  are  structured  into  a  format  ready  for 
DIDO  to  use  as  the  OCP  constraints.  Before  attempting  a  numerical  solution,  we  desire 
the  analytical  solar  sail  steering  law  to  help  verify  the  solution  as  described  below. 


C.  SOLAR  SAIL  CONTROL  LAW 

Following  the  guidelines  set  forth  by  the  minimum  principle  (Appendix  A),  we 
start  by  establishing  the  cost  function.  To  minimize  time,  we  choose  to  express  the  cost 
functional  from  the  previous  section  in  Lagrange  form. 

E  =  0,F  =1 


The  state  space  vector  and  its  derivativ  e  are  written  as 
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where  the  control  is  given  by  the  sail  cone  angle  a  (t) . 


Notice  that  0  is  uncoupled  from  the  other  equations,  therefore  it  will  be  ignored 
for  the  rest  of  the  sail  control  analysis.  Thus,  we  write  the  Hamiltonian  as 


H(x,u,X)  =  F  +  XT(  =  1  +  Xru  +  Xu 


where  X  is  the  costate  vector  and  the  solar  gravitational  parameter  has  been  normalized 
to  unity. 

The  Lagrangian  of  the  Hamiltonian  is  therefore 
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71  71 

where  u(/  is  now  the  covector  associated  with  the  path  constraint  —  —  <  a  <  — .. 


Applying  the  minimum  principle  yields 


dL 


^-=—X  ^-cos2asina  +  X 
da  "  r 


2(3 


— ■ — cos  oc  sin _  a  +  —cos  a 

2  2 

r  r 


+  H«  =0 


so  |i(/  may  be  written  as 

_  P 


(Xa  =-t-rcosa  (3A.ucosasina  +  2A.i,sin2a -X, cos2 a) 


For  the  given  control  constraints,  the  KKT  conditions  give 

a=“.Ha^0 

«=+|-^a  ^0 

—71  71  n 

—  <a<—  ,u  =0 
2  2  a 


— 71  71 

Limiting  the  analysis  to  the  interior  controls  where  — —  <a<  —  and  |i„=0,  we  can 


obtain  the  control  law. 
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( 3XU  cos  a  sin  a  +  2\,  sin 2  a  -  kv  cos 2  a )  =  0 

Dividing  the  equation  by  cos  2a  produces  a  quadratic  equation  with  tana  as  the 
independent  variable. 

2\  tan2 a  +  3ku  tana  -A,v  =  0 


(3) 


tana  = 


-3k  ±j9k  2  +8 A,2 

u  y  u  V 

4^ 


This  tangent  control  law  will  not  have  quadrant  ambiguity  since  we  limited  the 

control  angle  by  -  —  <  a  <  — . 

2  2 

D.  RESULTS 

Before  analyzing  results,  it  is  useful  to  correlate  certain  key  sail  angles  with 
physical  meaning.  As  mentioned  in  the  discussion  of  the  sail  model,  the  thrust  magnitude 
due  to  the  solar  radiation  pressure  (SRP)  on  the  sail  is  not  a  free  control  variable  as  it  is 
with  many  conventional  propulsion  systems.  The  acceleration  imparted  by  the  sail,  Ts/m, 
has  a  magnitude  that  is  dependent  on  radial  distance  and  the  control  angle,  a ,  according 

to  the  relationship  -  =  -^- cos 2  a  .  Notice  that  for  a  given  r,  maximum  radial  thrust 
m  r~ 

occurs  when  a  =  0  and  it’s  sail  area  is  totally  exposed  to  the  sun.  Alternatively,  the  sail  is 
effectively  “off’  and  the  trajectory  becomes  ballistic  when  a  =  ±71/2.  The  sail  can  never 
isolate  tangential  thrust  from  radial  thrust.  Radial  thrust  is  present  for  all  control  angles 
except  for  a  =  ±71/2  when  the  sail  edge  is  toward  the  sun.  To  assist  in  interpreting  the 
control  profile,  we  determine  the  control  angle  at  which  maximum  transverse  acceleration 
occurs.  The  transverse  acceleration  is  given  by: 

T.  .  Bu  2  7i  n 

a  =—  sina  cos~a  sina  ,  -  — <a<  — 

'  m  r2  2  2 

Setting  the  derivative  with  respect  to  a  to  zero  for  a  given  r  yields  the  control 

angle  providing  maximum  transverse  acceleration. 
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...  — L  =  -2cosa  sin2a  +  cos3a  =  0 

(4)  da 

a  =.615=  35.2° 

1.  Benchmark  Problem  Solutions  with  Circular  Coplanar  Orbits 
a.  Earth-Mars  Flyby 

The  purpose  of  solving  the  flyby  mission  serves  several  purposes.  First, 
the  solution  may  be  compared  to  other  known  solutions  to  verify  that  the  numerical 
analysis  using  DIDO  works  properly.  Second,  the  minimum  time  of  flight  serves  as  a 
lower  time  bound  for  future  coplanar  trajectory  optimization  problems.  Third, 
groundwork  is  established  to  model  a  bounded  initial  C3  from  Earth  allowing  a  li  mited 
“boost”  from  a  conventional  rocket  at  time  t  =  0. 

Given  an  initial  guess  of  the  initial  and  final  positions  and  velocities, 
DIDO  outputs  the  minimum  time  Mars  flyby  states  and  controls  (plotted  as  markers  in 
Figure  5).  For  a  minimum  time  flyby  path  between  circular  coplanar  orbits,  a  sail  with 
lightness  number  P=.  17  takes  0.45  years.  The  probe  sails  past  Mars  at  a  speedy  8.7  km/s 
(relative  to  Mars).  Sail  attitude  favors  a  large  local  transverse  acceleration  profile 
initially  to  rapidly  build  radial  acceleration,  and  then  gradually  exposes  more  sail  area  to 
the  sun  throughout  the  maneuver  until  the  sail  normal  is  parallel  with  the  sun’s  rays 
(Figure  6).  DIDO  costate  outputs  representing  X  -  [X( .,  \ ,  \ ,  k  ]7  are  plotted  in  Figure  7, 

and  may  be  used  with  the  tangent  steering  control  law  in  equation  (3)  to  generate  derived 
controls.  The  true  anomaly  costate,  \  ,  is  zero  as  expected  since  9  does  not  appear  in 
the  Hamiltonian.  A  comparison  of  costate -derived  controls  with  the  DIDO  output 
controls  appears  in  Figure  8. 
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Mars  min  time  flyby.  DIDO  vs.  propagated  states  and  controls.  16  nodes. 


Figure  5  States  and  Control  for  Mars  Flyby  Mission.  Markers  are  DIDO 
output  and  lines  represent  propagated  path. 
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Figure  6  Mars  Flyby  Trajectory  with  Sail  Profile  and  Normal. 


Costates  for  Mars  minimum  time  flyby  problem.  16  nodes. 


Figure  7  Mars  Flyby  Mission  Costates 
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Figure  8  Mars  Flyby  Cone  Angle  Control 

Thus  far,  the  flyby  model  has  restricted  the  sail  to  start  under  its  own 
power  without  the  aid  of  a  kick  motor.  Modeling  a  kick  motor  is  accomplished  by 
changing  the  initial  conditions  in  Table  2  for  the  u  and  v  velocity  components  to  be 
constrained  variables  according  to  the  following  relationship. 

0  <  IvJ  <  AV 

j  0  j  max 

where  V0  =  [n0,v0J7  —  V  .  Not  surprisingly,  the  optimal  solution  makes  use  of  whatever 
AVmax  is  permitted  to  intercept  Mars  in  the  quickest  time,  although  the  direction  of 
departure  depends  on  the  maximum  allowable  kick.  The  mission  is  accomplished  much 

faster  than  the  previous  sail-only  solution.  For  example,  when  AVmax  =  6  ------  the  optimal 

s 

solution  uses  all  the  initial  help  it  can  get  reducing  the  time  to  intercept  to  only  .22  years. 
We  find,  however,  that  the  same  behavior  does  not  apply  to  the  rendezvou  s  problem  with 
its  final  velocity  constraint. 

Feasibility  of  the  solution  is  demonstrated  in  Figure  9  depicting  the  path 
propagated  using  an  ODE  solver  given  the  initial  conditions  and  the  DIDO  -derived 
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control  profile.  Although  difficult  to  rigorously  prove  that  this  is  the  optimal  trajectory, 
we  can  at  least  show  that  necessary  conditions  for  optimality  are  satisfied  in  accordance 
with  the  discussion  in  the  section  on  Optimality.  The  Hamiltonian  is  constant  with 

respect  to  time  (Figure  10)  and  the  second  order  condition - ->  0  is  true  (Figure  11). 

da 

90  2 


Figure  9  Mars  Flyby  Mission  Propagation  (line)  with  DIDO  Output  (dots) 
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Figure  10  Hamiltonian  Profile  for  Mars  Flyby  Mission 


Figure  11  Second  Order  Condition  for  Mars  Flyby 
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b.  Rendezvous 

As  expected,  the  Mars  rendezvous  mission  takes  longer  to  accomplish 
than  the  flyby  mission,  requiring  a  full  1.11  years,  which  is  consistent  with  the  solution 
presented  in  ref  12,  (Figure  12).  The  spacecraft  final  velocity  is  constrained  to  match 
Mars’  orbital  velocity,  which  costs  the  sail  extra  transit  time  to  accomplish.  A  surprising 
feature  of  the  trajectory  is  that  the  probe  sweeps  an  arc  outside  of  the  circular  orbit  of 
Mars  (Figure  13).  This  maneuver  is  observed  in  more  complex  minimum  time  solar  sail 
missions  as  well  since  it  appears  to  make  the  best  use  of  radial  and  transverse  thrust. 
There  is  a  point  on  the  outbound  path  (almost  a  quarter  of  the  way  into  the  mission  time) 
when  the  control  angle  rotates  from  a  small  angle,  where  much  of  the  sail  area  is  exposed 
to  the  sun,  to  a  more  edge-on  aspect.  The  sail  attitude  gradually  rotates  after  this  to  favor 
more  transverse  acceleration  as  it  sweeps  past  Mars’  orbit  radius  and  approaches  Mars. 
By  the  time  the  sail  has  reached  Mars,  the  control  angle  is  5°  greater  than  its  maximum 
transverse  velocity  setting. 


EM  Rendezvous  with  50  nodes 


Figure  12  DIDO  State  and  Control  Output  (Markers)  and  Propagated  Path 

(line)  for  EM  rendezvous. 
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EM  rendezvous  with  50  nodes 


Figure  13  Sail  Trajectory  and  Profile  for  Earth-Mars  Rendezvous 


Interpreting  the  costate  history  provides  a  more  complete  physical  analysis 
of  the  optimal  trajectory.  Each  costate  in  Figure  14  represents  a  Lagrange  multiplier  that 
signifies  instantaneous  sensitivities  of  the  cost  function  to  inst  antaneous  variations  in  the 
corresponding  state.  Once  again,  the  costate  corresponding  to  0  is  zero  since  this  state  is 
completely  uncoupled  from  the  other  states  in  the  dynamic  equations  and  never  appears 
in  the  Hamiltonian.  Several  observations  may  be  made  regarding  the  critical  point  in  the 
first  quarter  of  the  trajectory  when  the  control  angle  makes  a  radical  change.  Recalling 
the  optimal  sail  steering  law  in  equation  (3)  it  is  seen  that  when  \  =0,  then  the 


maximum  transverse  velocity  profile  occurs.  This  value  would  occur  when  tana  = 


£ 

2 


(a  =35.2°)-  Also,  when  \  — >0  then  a  —>90°.  Observing  the  costates  for  the  Earth- 


Mars  (EM)  rendezvous,  we  note  that  \  crosses  zero  at  t  =  1.225  years  and  soon 
afterwards  Xv  reaches  its  closest  point  to  zero.  This  is  consistent  with  a  rapid  change  of 
sail  attitude  at  that  time  during  the  mission  and  corresponds  to  when  the  sail  reaches 
maximum  radial  velocity. 
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Lagrange  multipliers  along  with  state  outputs  are  fed  into  the  derived 
tangent  steering  law  (equation  (3))  to  obtain  a  control  profile.  The  resulting  control 
profile  is  shown  in  Figure  15  and  is  plotted  with  the  DIDO  control  output  for  comparison. 


Costates  for  EM  rendezvous  mission.  50  nodes. 


Time  (years) 

Figure  14  Costates  for  Earth- Mars  Rendezvous  Mission 
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Control  Profile 


Figure  15  Comparison  of  DIDO  Controls  and  Tangent  Steering  Control  Law  for 

EM  Rendezvous. 


Propagating  the  spacecraft  with  the  given  control  history  produces  the  path 
shown  in  Figure  17  with  DIDO  states  at  node  points  shown  for  comparison.  Testing  the 


d2H 


second  order  necessary  conditions  for  optimality  we  observe  in  Figure  16  that - -  >  0 

dor 


for  all  time  and  also  the  Flamiltonian  is  fairly  constant  (Figure  18). 
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oH 

Figure  16  Second  Order  Necessary  Condition  — -  >  0 

3a  “ 


90  2 


Figure  17  Propagated  Path  (line)  with  DIDO  State  Output  (dots)  for  EM 

Rendezvous 
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50  nodes. 


Time  (years) 

Figure  18  Hamiltonian  for  Minimum  Time  EM  rendezvous. 


What  happens  when  we  reverse  the  problem  such  that  we  seek  the 
minimum  time  from  Mars  (VJmars  =  0)  to  Earth  rendezvous?  It  turns  out  that  swapping 
initial  conditions  with  final  conditions  generates  a  time -reversed  state  history  and  a 
control  profile  that  is  inverted  and  reversed.  The  experiment  may  be  earned  out  within 
the  DIDO  structure  by  simply  reversing  the  event  conditions. 
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Event 

ME  Rendezvous 
[lb,ub] 

ro 

[ a  ,a  ] 

1  m 7  m  J 

00 

free 

U0 

[0,0] 

vo 

1  1 

rf 

[1,1] 

0/ 

free 

uf 

[0,0] 

vf 

[1,1] 

Table  3  Mars-Earth  Rendezvous  Event  Conditions. 

The  state  history  of  the  optimal  Mars-Earth  rendezvous  is  a  mirror  image 
of  the  Earth-Mars  rendezvous  (Figure  19  and  Figure  20).  Time  of  flight,  At,  and  the 
change  in  angular  displacement,  A0,  are  precisely  the  same  for  both  profiles. 
Furthermore,  the  control  profile  is  inverted  and  time  reversed.  This  result  leads  to  an 
even  more  interesting  question;  what  happens  when  we  desire  the  sail  to  launch  from 
Earth,  rendezvous  with  Mars,  then  immediately  start  its  trek  back  home  to  rendezvous 
with  Earth?  This  benchmark  mission  is  the  subject  of  the  next  section. 
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Solar  sail  trajectory  with  sail  normal  vector  history 


Figure  19  Mars  to  Earth  Rendezvous.  Shows  reverse  trajectory  of  EM 

rendezvous  mission 


normalized  time  units 


Figure  20  ME  Rendezvous  States  and  Controls 
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c. 


Double  Rendezvous 


As  a  next  step  toward  modeling  a  cycler  orbit,  we  explore  the  double 
rendezvous.  In  this  mission  we  seek  the  minimum  time  that  it  takes  to  make  an  Earth  - 
Mars-Earth  (EME)  round  trip  with  a  zero  relative  velocity  at  both  planets  for  a  given  sail 
lightness  number  p.  To  study  the  behavior  of  the  optimal  trajectory,  one  assumes  that 
Earth  or  Mars  gravity  will  not  significantly  influence  the  craft.  Low  relative  velocity  is 
desirable  while  encountering  the  target  planets  since  it  requires  less  energy  for  a  greeting 
taxi  craft  to  intercept  and  dock  with  the  passing  spacecraft. 

For  the  preliminary  analysis,  the  initial,  intermediate  and  final  orbits  are 
modeled  as  circular  and  coplanar.  Earth  phasing  for  the  final  Earth  rendezvous  event  is 
accomplished  by  coupling  the  sail  and  Earth  final  positions  using  the  Earth’s  mean 
motion  and  the  time  of  flight. 

The  solution  to  the  minimum  time  EM  rendezvous  problem  provides  a 
lower  bound  for  the  time  that  it  must  take  to  reach  Mars  during  the  double  rendezvous 
problem.  Since  the  fastest  possible  time  to  rendezvous  with  Mars  with  the  given  solar 
sail  (P=.17)  was  determined  to  be  1.11  years  from  the  previous  section,  this  serves  as  the 
minimum  time  to  encounter  Mars  halfway  through  the  double  rendezvous  mission. 
Remembering  that  the  reverse  Mars-Earth  (ME)  rendezvous  took  exactly  the  same 
amount  of  time,  we  now  have  bounded  the  time  to  complete  the  whole  EME  double 
rendezvous.  Setting  the  intermediate  destination  at  the  Mars  orbit  and  the  final 
destination  with  Earth  and  constraining  the  event  velocities  to  match  those  of  the 
respective  planets  we  are  in  a  position  to  solve  the  optimal  EME  double  rendezvous. 
Shown  in  Table  4  is  a  summary  of  the  event  constraints  where  superscripts  indicate 
before  (-)  or  after  (+)  the  intermediate  knot  representing  the  Mars  encounter.  There  is  a 
single  non-linear  event  constraint  that  is  responsible  to  ensure  proper  phasing  (i.e.  Earth 
makes  more  orbits  around  the  sun  than  the  sail  does,  but  still  meets  with  it  in  the  end). 
The  output  states  and  controls  are  shown  Figure  21  with  the  path  shown  in  Figure  22. 
Essentially  it  turns  out  to  be  the  EM  minimum  time  rendezvous  solution  patched  together 
with  the  ME  rendezvous  solution  with  slightly  longer  transit  times.  The  most  obvious 
features  of  the  profile  are  the  state  symmetries  and  control  antisymmetry  about  the  mid  - 
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maneuver  point.  Event  conditions  at  the  interior  knot  required  only  that  the  states  be 
continuous  and  equal  to  Mars’  state. 
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Table  4  Event  Conditions  for  EME  Double  Rendezvous 


For  the  given  sail,  the  total  time  to  complete  the  trajectory  is  2.41  years 
with  both  legs  of  the  journey  taking  1.205  years.  The  reason  that  the  EME  problem  takes 
longer  than  the  patched  EM -ME  problem  is  that  the  final  Earth  rendezvous  event  must  be 
phased  with  Earth’s  position  at  the  final  time.  The  individual  EM  and  ME  rendezvous 
solutions  only  served  to  provide  a  lower  time  bound  for  the  EME  pr  oblem,  so  imposing 
the  constraint  that  the  sail  final  position  is  phased  with  Earth  is  about  10%  more  than  the 
lower  bound. 
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EME  rendezvous  -  min  time 


Figure  21  States  and  Controls  for  EME  Double  Rendezvous  in  Minimum  Time 

(P=-17) 


Figure  22  EME  Double  Rendezvous  Trajectory  and  Sail  Profile  (P=.17) 
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This  round  trip  solution  only  takes  into  account  a  quick  trip  out  to  graze 
the  Mars  orbit  and  then  a  return  trip  home.  It  may  be  desirable  to  have  the  probe  linger 
with  Mars  for  a  span  of  time  allowing  for  some  sort  of  mission  while  traveling  with  the 
planet  then  returning  back  to  Earth.  The  next  question  one  may  ask  is  the  following. 
Given  a  sail  of  certain  performance  and  Mars  stay  time,  what  is  the  minim  um  transit  time 
EME  trajectory  (not  including  the  Mars  stay  time)?  For  a  given  Mars  stay  time,  some  of 
the  transit  time  to  and  from  Mars  is  expended  in  performing  a  phasing  maneuver  to  meet 
Earth  at  the  final  time.  This  characteristic  phasing  maneuver  is  exhibited  in  all  double 
rendezvous  problems  with  varying  stay  times  but  one  -  the  one  where  the  stay  time  spent 
with  Mars  happens  to  be  just  the  right  time  to  allow  a  return  trip  to  Earth  with  no  phasing 
maneuver  required  at  all.  As  mentioned,  this  minimum  possible  return  trip  time  is  1.1 
years  for  sail  lightness  number  P=.  17.  The  following  analysis  shows  that  given  the 
optimal  time  and  angular  displacement  required  for  the  EM  rendezvous  (and  thus  the  ME 
rendezvous),  we  can  deduce  the  Mars  stay  time  that  minimizes  transit  time.  The  optimal 
transit  time  and  true  anomaly  traversed  by  the  sail  for  an  EM  rendezvous  (and  ME 
rendezvous)  are  given  by  the  following. 

t*  =  1.116  years  0*  =  4.337  radians 

The  following  equations  c  apture  the  motion  of  the  Earth  relative  to  that  of 
the  cycler  craft  where  tm=  Mars  stay  time,  A0  =  change  in  Earth  angular  position  and 
ne  and  nm  represent  Earth  and  Mars  mean  motion  respectively. 

(5)  A0e  =  ne(2f  +  tm) 

Equation  (5)  indicates  that  Earth  changes  position  by  A0f  during  the  same 
time  that  the  sail  transits  outbound  in  t  years,  stays  with  Mars  for  tm  years,  then  returns 

home  in  t  years.  Recognizing  the  n  we  can  rewrite  equation  (5)  as 

(6)  A0  -t  =2 1* 

v  z  e  m 
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During  the  round  trip  mission  time  the  spacecraft  traverses  0  *  radians  on 
the  outbound  leg,  nmtm  radians  along  with  Mars,  and  0  *  on  the  inbound  leg. 

(7)  A0  =20*  +  «  t 

v  '  mm 


All  the  while,  Earth  traverses  the  same  span  plus  an  extra  N  revolutions. 

(8)  A0f  =  A0  +27tiV£  (Ne  =[0,1,2...]) 


Constraining  the  solution  to  include  only  the  first  Earth  pass  on  the  return  (N=l),  we 
substitute  equation  (7)  into  equation  (8)  and  rearrange  to  obtain 

(9)  A0c  —nmtm=  20*+  2k 


Recognizing  that  the  scaled  mean  motions  of  Earth  and  Mars  are  given  by 


n,  =1  and  n,„  =  ./ — -  =  .534  respectively,  we  may  solve  for  A0f 


using  equations  (6)  and  (9). 
days  and  an  Earth  span  of 
mission  time. 


This  results  in  a  Mars  stay  time  of 
A0f  =  16.03  radians  corresponding 


and  tm  simultaneously 

tm  =7.0114TUs  =  117 
to  2.55  years  of  total 
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Figure  23  EME  Double  Rendezvous  with  Stay  Time  at  Mars 

When  using  the  outbound  and  inbound  transit  times  as  a  cost  function, 
DIDO  produces  the  output  shown  in  Figure  23.  The  path  mimics  a  patched  minimum 
time  EM-ME  solution  where  the  optimal  Mars  stay  time  is  123  days  for  the  given  sail 
(5%  difference  from  the  estimated  value  of  1 17  days).  Any  more  or  less  time  spent  at  the 
planet  will  either  require  a  time-wasting  phasing  maneuver  to  allow  Earth  to  catch  up  for 
a  rendezvous,  or  a  maneuver  to  catch  up  to  the  speeding  target  Earth.  A  plot  showing  the 
impact  of  various  Mars  stay  times  on  mission  transit  times  appears  in  Figure  24.  Notice 
that  it  is  safer  to  design  a  shorter-than-optimal  Mars  stay  time  into  the  mission.  In  the 
event  of  a  schedule  slide,  the  return  transit  time  looks  increasingly  grim  beyond  the 
optimal  point  because  the  sail  has  a  lot  of  space  to  cover  as  it  attempts  to  catch  up  to 
Earth. 
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Figure  24  EME  Double  Rendezvous  Transit  times  for  Various  Mars  Stay  Times. 

((3  =.17) 

The  analysis  thus  far  has  only  dealt  with  sails  of  lightness  number  (3=.  17. 
An  interesting  question  that  may  be  asked  now  is  the  following.  Given  a  desir  ed  time  of 
round  trip  flight  with  no  stay  time  at  Mars,  what  is  the  minimum  size  solar  sail  required? 
This  question  will  be  addressed  later. 

2.  Benchmark  Problem  Solutions  with  an  Elliptic  Coplanar  Orbit 

Running  the  same  battery  of  problems  from  the  previous  section  using  a  higher 
fidelity  model  provides  enormous  insights  into  the  characteristics  of  an  optimal 
trajectory.  Modeling  trajectories  between  circular  Earth  and  coplanar  elliptic  Mars  orbits 
generates  results  that  may  be  compared  with  the  circula  r  orbits  model  to  reveal  the  “knees 
in  the  curve”  that  the  optimization  process  seeks  in  an  effort  to  reduce  the  cost  function 
just  a  little  more.  Knowing  these  characteristics  of  optimal  solar  sail  trajectories  will 
assist  in  understanding  the  behavior  of  optimal  solar  sail  cyclers. 

Because  the  final  target  manifold  represents  Mars’  elliptical  orbit,  the  new  events 
model  must  relate  the  final  Mars  radial  position,  rmf ,  and  velocity,  \m ,  with  the  final  sail 
state,  x  . 
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Constraining  the  initial  and  final  sail  angular  positions  to  be  coincident  with  the 
initial  Earth  and  final  Mars  angular  positions  respectively,  we  obtain  the  following 
relationship  for  coplanar  orbits. 

9„,  =e/-(®„+*U 

where  0;  is  the  final  sail  position,  com  is  Mars’  argument  of  periapsis  and  £lm  is  the  right 
ascension  of  the  ascending  nodes. 

In  a  perifocal  system,  Mars  final  polar  coordinates  rmf  and  0  are  related  by 


(10) 


'  mf 


l  +  e  cos0 


mf 


where  pm  in  the  numerator  is  the  “parameter”  or  semi-latus  rectum  of  the  Mars  orbit. 
The  final  event  condition  for  radial  position  of  the  sail  is  simply 


(ID 


r.  =  r. 

f  if 


For  the  rendezvous  mission,  we  need  to  target  the  final  Mars  velocity  as  well.  The  speed 
of  Mars  at  the  final  time  is  expressed  as 


r  . 
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[2  1  | 
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a  1 
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\  mf 

m  1 

^ «f  I1 

where  am  is  the  semi-major  axis  of  Mars.  Components  of  the  planetary  velocity  vector 
are  defined  in  terms  of  the  flight  path  angle  (3m  by 


cos  Pm/  =  — and  sin  \i,rf  =-^—em  sin  (0m/ ) 

Cl/  mf  m  nf 

where  hm  represents  the  magnitude  of  the  Mars  angular  momentum  vector. 

The  velocity  vector  resolved  in  the  local  vertical  local  horizontal  frame  (LVLH) 
produces 

Unf  =K,/SinP„/ 
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v  =  V  .  cos  B  , 

mj  try  *  nf 

Representing  final  Mars  and  sail  velocities  by  V  =  [umf  ,vmf  ]T  and  V  =  [uf  ,vf  f 
respectively,  we  obtain  an  additional  constraint  for  the  elliptic  rendezvous  mission 
(12)  V,  =  V  , 

It  now  remains  to  run  the  same  battery  of  battery  of  benchmark  problems  with  the 
elliptic  Mars  event  model. 


a.  Earth-Mars  Flyby 

Modeling  Mars’  orbit  as  an  ellipse  using  eccentricity  em  =  .0935  provides 
target  radial  distances  that  vary  with  the  planets  true  anomaly  (equation  (10)).  Earth’s 
orbit  eccentricity  is  only  .0167,  so  the  circular  Earth  orbit  assumption  is  a  good  one 
(elliptic  Earth  orbit  is  considered  in  the  3-D  model  in  the  next  chapter).  Not  surprisingly 
the  time  optimal  mission  selects  a  path  that  drives  the  sail  from  Earth  orbit  to  Mars 
periapsis  at  high  speed  (Figure  25).  This  optimal  path  exploits  Mars’  eccentric  orbit  and 
presents  the  shortest  distance  from  Earth  to  Mars,  and  thus  the  shortest  transit  time 
required  since  the  final  velocity  is  unconstrained.  The  time  to  intercept  Mars  using  the 
standard  sail  (()  =.17)  is  only  137  days,  about  27  days  faster  than  the  corresponding 
circular  orbits  model  (Figure  26). 
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1 6  nodes. 


Figure  25  Minimum  Time  EM  flyby  with  Mars  Elliptical  Orbit.  Final  position 

is  at  Mars  periapsis. 


EM  minimum  time  flyby  with  elliptic  coplanar  orbits.  16  nodes. 


Figure  26  EM  Flyby  with  Mars  Elliptic  Orbit 
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b. 


Rendezvous 


To  achieve  the  rendezvous  mission  in  minimum  time  the  path  must 
account  for  target  velocity  as  well  as  distance,  both  of  which  vary  with  respect  to  target 
planet  true  anomaly.  Initial  and  final  true  anomalies,  thus  launch  windows  and  arrival 
dates,  are  optimization  parameters  that  seek  to  reduce  the  cost  function  by  finding  the 
optimal  points  on  the  target  manifolds  to  bound  the  path.  It  is  fascinating  to  observe  that 
the  optimal  rendezvous  with  Mars  (elliptic  orbit)  is  actually  faster  than  the  simplified 
circular  coplanar  rendezvous  (0.977  years  vs.  1.11  years).  Comparison  of  this  state 
history  with  the  circular  coplanar  counterpart  gives  insights  as  to  why  this  occurs  (  Figure 
27).  In  the  circular  coplanar  model,  any  mission  start  date  is  as  good  as  any  other 
therefore  the  boundary  condition  0O  =  0  was  a  valid  restriction  where  the  Mars  lead  angle 
was  determined  in  the  numerical  solution  to  the  OCP.  Because  in  the  new  events  model 
we  allowed  the  optimal  0()  to  be  determined  in  the  OCP  solution,  a  launch  window  was 

chosen  such  that  Mars  would  be  near  the  slowest  point  in  its  elliptic  path  at  the  time  of 
rendezvous.  Mars  matches  the  sail  velocity  at  the  final  event  as  it  slows  down  on  its 
approach  to  aphelion.  The  optimal  rendezvous  trajectory  takes  advantage  of  Mars’  slow 
velocity  as  the  sail  approaches  the  planet  at  a  radial  distance  that  is  95%  of  the  aphelion 
distance.  Figure  28  shows  how  the  sail  meets  the  Mars  position  on  its  orbit. 

Another  interesting  departure  of  the  elliptic  orbit  solution  from  the  circular 
solution  is  that  the  spacecraft  does  not  follow  a  trajectory  that  sweeps  out  to  a  maximum 
radius  and  then  returns  inward  to  meet  Mars  at  the  required  velocity  (Figure  27).  Mars’ 
orbital  path  however  does  cross  inside  the  spacecraft  path  as  Figure  28  readily  reveals. 
The  optimal  3-D  trajectory  takes  the  fastest  path  from  the  optimal  location  on  Earth’s 
orbit  to  intercept  Mars  without  having  to  undergo  a  negative  radial  velocity  during  the 
whole  maneuver,  i.e.  u  >  0  for  all  time  (Figure  29). 
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Optimal  state  paths  for  elliptic  and  circular  orbits 


Figure  27  Comparison  of  State  Profiles  Using  the  Circular  (markered  lines)  and 
Elliptic  (thick  lines)  Orbits  Models. 

40  nodes 


Figure  28  EM  Rendezvous  with  Mars  Elliptic  Orbit 
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40  nodes 


c.  Double  Rendezvous 

Once  again,  aphelion  provides  the  optimal  Mars  rendezvous  point  ( Figure 
30  and  Figure  1).  The  total  transit  time  is  only  2.3  years  in  contrast  to  the  24  years  it 
took  to  reach  the  circular  Mars  orbit.  In  contrast  to  the  double  rendezvous  orbits  (see 
Figure  21),  the  radial  velocity,  u ,  does  not  change  as  much  in  the  vicinity  of  the  interior 
knot  thus  “flattening”  the  radial  distance  profile  as  the  spacecraft  sweeps  past  Mars. 
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EME  double  rendezvous  between  elliptic  orbits.  50  nodes. 


Figure  30  States  and  Control  for  EME  Double  Rendevous  with  Elliptic  Mars 

Orbit. 


50  nodes. 


Figure  31  EME  Double  Rendezvous  Trajectory  and  Sail  Profile  with  Elliptic 

Mars  Orbit. 
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3. 


Earth- Mars  Synodic  Cycler  Solution 


Having  established  a  collection  of  solutions  for  basic  solar  sail  trajectories,  we  are 
in  a  position  to  take  the  next  step  toward  a  cycler  model.  There  are  two  key  differences 
between  the  double  rendezvous  mission  of  the  previous  section  and  the  solar  sail  cycler 
mission.  With  cyclers,  end  conditions  are  equal  (i.e.  cyclic)  by  definition  to  ensure 
repeatability.  Additionally,  gravity  assists  at  target  planets  will  be  modeled.  Planetary 
swingbys  offer  enhanced  cycler  performance  since  the  turn  angles  are  optimized  to  create 
a  path  that  achieves  the  minimum  cost.  Since  the  time  to  complete  a  cycle  is 
predetermined  by  the  synodic  period  (Figure  32),  a  minimum  time  synodic  cycler 
problem  has  no  meaning  for  spacecraft  shuttling  between  circular  coplanar  orbits.  This  is 
not  so  for  the  cycler  between  non -circular  coplanar  orbits.  A  new  cost  function  is  desired 
that  is  not  burdened  with  minimizing  fuel  or  time.  Because  an  attractive  feature  of 
cycling  trajectories  is  a  slow  swingby  velocity  at  Earth  and  Mars,  these  velocities  formed 
the  cost  function. 


a.  Events  Model 

The  key  to  modeling  a  repeating  cyclic  trajectory  is  to  constrain  the  initial 
state  of  the  spacecraft  to  equal  the  final  state  of  the  spacecraft.  Since  the  Earth  and  Mars 
orbits  are  approximated  as  circular  and  coplanar,  the  problem  is  simplified  in  that  the 
initial  relative  angular  position  between  Earth  and  Mars  (lead  angle  £  )  may  be  used  to 
constrain  final  angular  position  for  a  single  cycle  since  planetary  distances  and  velocities 
are  independent  of  the  their  inertial  angular  positions  (Figure  32).  Lead  angle  is 
determined  from  the  event  s  tates  as  follows. 


=  0  n  +  n 

ml)  n 


( tf  -t0) 


—  ne^tf  ~ 

where 

0  n  =0  —  n  (t  -tn) 

ml)  i  m  k  i  O' 
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Mars  and  Earth  angular  (0m  and  0  )  are  obtained  using  their  respective 


mean  motions  (n  and  n  ). 

v  m  e ' 


The  Earth-Mars  lead  angle  at  initial  and  final  times  are 


—  9»e 

C/  =®f,f  ~®ef 


Initial  and  final  relative  angular  positions  are  constrained  by 

(13)  cos(C/-C0)  =  l 

where  C,0  =  C>f -2kN^  ,  N r?  =0,1,2,... 

Numerically,  equation  (13)  is  preferable  over  the  above  equation  since  the 
cosine  function  will  permit  multiple  revolutions  without  introducing  integer  variables. 

The  cyclic  end  condition  for  the  spacecraft  radial  distance  is  expressed 

simply  as 

(14)  r0=rf 

Initial  and  final  velocities  are  also  constrained  to  be  cyclic  and  will  be 
addressed  later. 

Given  the  condition  in  equation  (13),  the  final  time,  tf ,  depends  only  on 

the  relative  mean  motions  of  the  planets  given  by  the  synodic  period, 

2n 

x<  =i - 

\n  -n 


For  Earth  and  Mars,  the  synodic  period  is  2.135  years'1.  Although  not 
explicitly  constrained  to  the  synodic  period  in  the  numerical  analysis,  the  resulting  final 
time  must  equal  this  synodic  period  in  order  for  the  cycler  to  be  periodic.  Having 
established  cyclic  end  conditions,  we  now  turn  to  setting  up  gravity  assist  conditions  for 
planetary  encounter  events. 
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Figure  32  Two-Dimensional  Earth-Mars  Cycler  Geometry  with  Circular 

Coplanar  Planetary  Orbits 

Force  imparted  by  the  solar  sail  shapes  the  path  between  event  manifolds, 
but  gravity  assist  maneuvers  at  the  event  manifolds  drive  the  form  of  the  whole  cycler 
trajectory.  As  with  the  conventional  Aldrin  cycler,  gravity  assists  are  implemented  in  the 
solar  sail  cycler  to  shape  the  trajectory  at  these  planetary  encounter  events  to  improve  the 
revisit  times.  To  model  the  swingby  events,  state  discontinuities  are  employed  at  the 
planet  to  change  the  velocity  direction  and  magnitude  of  the  spacecraft  in  the  heliocentric 
frame  (sometimes  called  the  “zero  sphere  patched  conic”14  or  matched  asymptotes 
model).  It  is  assumed  that  the  interaction  time  with  the  subject  planet  is  negligibly  small 
in  comparison  to  the  total  cycle  time.  The  velocity  of  the  spacecraft  with  respect  to  the 
planet  changes  direction  such  that  Vtp+  =  V  “+AV,  where  Vip  =  [ u,v]r  —  V  and  Yp  is 

the  velocity  of  the  subject  planet  relative  to  the  sun.  The  position  states,  r  and  9  are 
constrained  to  be  continuous  at  both  Earth  (equation  (14))  and  Mars  encounter  events 
given  below. 

r=  r+ 

(15) 

e~=e+ 
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The  AV  s  due  to  the  swingbys  are  optimally  chosen  in  the  OCP  solution, 
but  they  must  be  properly  constrained.  Constraints  on  the  velocity  changes  are  most 
easily  imposed  using  the  planet  frames  where  inbound  and  outbound  velocity  magnitudes 
are  equal  just  before  and  after  a  planet  encounter  forming  the  event  condition, 

(16)  |V  |  =  |V„+ 1  =  K.  (planet  frame) 

where  VC  is  the  hyperbolic  excess  speed.  The  velocity  direction  change  is  expressed  in 

the  planetary  frame  using  the  turn  angle,  8  ,  which  exists  in  the  region  shown  in  . 
Velocities  before  and  after  a  swingby  event  in  the  planet  frame  are  coupled  by  the  cosine 
of  the  turn  angle, 

V  -.y  +  =v  2cosS 

sp  sp  oo 

V  *v  + 

(17)  cos  5  =  sPy  , sp  (planet  frame) 

The  spacecraft  experiences  a  direction  change  during  the  interaction  that  is 
restricted  by  the  hyperbolic  excess  speed  and  the  permissible  periapsis  pass  distance  from 
the  center  of  the  subject  planet,  rp .  This  restriction  is  expressed  by  the  following 

relationship  in  the  planet  frame  where  VC  is  scaled  by  the  circular  orbit  speed  at  the 
surface  of  the  planet  and  rp  is  scaled  by  the  radius  of  the  planet  (ref  15  p.  24). 

(18)  sin (8  /2)  = — - - ,  r  .  <r  <r 

v  '  i  V~r  +1  pnun  p  pmax 

“  p 

Although  equations  (17)  and  (18)  themselves  are  not  event  conditions, 
they  limit  the  achievable  change  in  velocity  direction,  AV ,  due  to  the  swingbys  at  Mars 
and  Earth.  Design  limitations  include  a  selected  rpmm  that  is  well  above  the  atmosphere 
where  drag  effects  are  negligible  and  an  r  such  that  the  encounter  occurs  close  enough 
to  the  planet  to  execute  a  desired  task.  For  a  given  VC ,  the  turn  angle  is  maximum  when 
r  =  rpmin  and  minimum  when  rp  =  rpmax .  Substituting  these  values  into  equation  (18)  and 
solving  for  8  as  a  function  of  VC  provides  an  expression  for  8maxand  8min  respectively. 
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To  include  the  case  of  maximum  deflection  in  the  opposite  direction  it  is  necessary  to 
place  the  lower  bound  of  acceptable  turn  angles  at  -Sminand  -5max  (Figure  33).  With 
equation  (17)  characterizing  the  instantaneous  change  in  path  direction  and  equation  (18) 
providing  the  limits,  the  boundary  conditions  are  expressed  at  Earth  and  Mars  event 
manifolds  as 

(19)  cos(8max)  <cos8  <cos(8min) 


Figure  33  Gravity  Assist  Geometry. 


Velocity  constraints  at  Mars  and  Earth  have  identical  form,  however  the 
initial  velocity  magnitude  at  Earth  has  an  additional  limitation  based  on  available 
departure  rocket  capability.  Because  the  sail’s  journey  starts  at  Earth,  the  initial 
conditions  are  bounded  by  maximum  C3  available.  Presumably  an  impulse  rocket  is  used 
to  start  the  solar  sail  craft  on  its  cycler  trajectory,  so  the  initial  velocity  relative  to  Earth, 
V0  \earth ,  is  restricted.  The  magnitude  of  V0  \earth  is  limited  by  a  maximum  allowable 

velocity  change  at  t  =  0  provided  by  a  kick  motor,  which  provides  another  boundary 
condition 


(20) 


0  ^  |Vo  Lnk  |  ^  AVmax  (Earth  frame) 
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The  direction  of  V0 1  k  is  driven  by  the  optimal  control  problem  and  is 
limited  only  by  the  allowable  turn  angle  at  each  Earth  swingby. 

Finally,  phasing  the  spacecraft  with  Earth  at  the  end  of  a  cycle  was 
considered.  Earth  encounter  events  were  constrained  to  ensure  that  the  sail  trajectory 
intersected  Earth’s  orbital  path  at  precisely  the  time  that  the  planet  is  at  that  same 
location.  The  circular  orbit  assumption  is  particularly  useful  for  ensuring  proper  phasing 
of  events  since  the  angular  position  of  a  planet  is  a  linear  function  of  time.  The  final 
angular  position  is  given  by 

(21)  0y  =0eO  +ntf  -2nNs, 

where  Ns  is  an  integer  number  of  Earth  orbits,  ne  is  Earth’s  mean  motion,  and  0eO  is  the 
angular  position  of  Earth  at  t  =  0  (when  we  the  Earth  position  is  coincident  with  the  sail, 
0(n  =0O).  Having  established  the  sail,  dynamic  and  events  models,  we  set  up  the  optimal 
control  problem. 

b.  Solar  Sail  Cycler  Problem  Formulation 

The  solar  sail  cycler  optimal  control  problem  is  constructed  using 
weighted  spacecraft  Ws  at  the  Mars  and  Earth  encounters  as  the  cost  function.  The 
optimal  path  is  subject  to  two-dimensional  equations  of  motion,  cyclic  end  conditions, 
and  planetary  gravity  assist  constrain  ts.  The  cost  function  uses  a  parameter  y  to  weight 
the  Vo„s  while  the  initial  and  final  state  and  control  conditions  are  constrained  to  be  equal 
to  ensure  repeatability  of  the  trajectory.  We  can  write  the  optimal  control  problem  as  the 
following. 


Minimize  the  cost 


(22) 


7(x,,y)=yVJ 


+d-Y)V  u 


Subject  to  dynamic  constraints  from  the  equations  of  motion  (equation  (2)) 
and  event  constraints  that  model  repeatability  (equations  (13)  and  (14)),  continuity 
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(equation  (15)),  swingby  effects  (equations  (16)  and  (19)),  launch  limitations  (equation 
(20))  and  phasing  (equation  (21)). 

In  addition  to  the  bounds  on  controls  (equation  (1)),  there  are  bounds  on 
states.  States  are  bounded  only  to  avoid  singularities  in  the  dynamics. 

For  this  paper,  a  single  Earth-Mars  cycle  was  modeled  where  the 
following  initial,  intermediate  and  final  conditions  in  astronomical  units  define  the  target 
manifolds 

ro  =rf  =  1 
r  =1.524 

The  parameters  that  bound  the  path  deflections  at  the  event  manifolds  are 
the  maximum  Q  available  at  the  initial  Earth  orbit  departure,  AVmax  in  canonical  units, 
and  the  minimum  and  maximum  periapsis  pass  distances  at  Earth  and  Mars. 

AVmax=.2  (~6 km/s) 

r  =  1 .06  Mars  radii 

mmin 

r  —  unbounded 

m  max 

rmin  =1.16  Earth  radii 

r  =10  Earth  radii 

<?max 

c.  Synodic  Cycler  Results  and  Analysis 

Shown  in  Figure  34  is  the  state  and  control  angle  output  from  DIDO  for  a 
single  cycle  of  the  synodic  cycler  with  y=l.  As  expected,  the  time  required  to  complete  a 
cycle  under  this  set  of  constraints  was  2.135  years,  the  Earth-Mars  synodic  period.  A 
quick  glance  at  the  state  history  reveals  that  the  spacecraft  sails  from  1  AU  out  to  Mars’ 
orbit  at  1.524  AU  and  then  returns  to  Earth.  A  discontinuity  in  both  velocity  states  occurs 
at  1.524  AU  and  1  AU  representing  Mars  and  Earth  swingbys  respectively.  The  optimal 
path  makes  use  of  large  gravity  assist  maneuvers  (in  the  planet  frame)  during  planet 
encounters  owing  to  slow  Vo„s.  The  spacecraft  initially  accelerates  in  the  radial  direction 
while  decelerating  in  the  transverse  direction.  At  the  appropriate  time,  the  sail  rotates  to 
an  attitude  that  favors  more  and  more  positive  transverse  acceleration  to  intercept  Mars 

with  the  lowest  Y*,  to  minimize  the  cost  function  while  setting  up  for  the  swingby  event 
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initiating  the  return  trip.  Following  Mars  swingby,  the  spacecraft  sweeps  out  to  nearly  2 
AU  to  ensure  proper  phasing  for  Earth  intercept.  Sail  attitude  gradually  reaches  a 
maximum  negative  transverse  acceleration  profile  (a  =  -35°;  see  equation  (4)),  then 
“shuts  off’  and  follows  a  ballistic  path  as  it  presents  an  edge  aspect  to  the  sun.  A  plot  of 
the  sail  trajectory  with  Earth  1 ,  Mars  2  and  Earth  3  encounters  is  shown  in  Figure  35.  A 
similar  gravity  assist  is  accomplished  at  the  Earth  encounter  and,  since  cyclic  end 
conditions  were  imposed,  the  same  control  profile  will  reproduce  the  trajectory 
repeatedly.  Owing  to  these  constrained  end  conditions,  the  initial  Earth  departure 
hyperbolic  excess  velocity  only  required  4.3  km/s  -  not  the  maximum  allowable  limit  of 
6  km/s. 


Synodic  Cycler  with  1 15  total  nodes 


Figure  34  DIDO  States  (markers)  and  Control  with  Propagated  Path  (line 
through  markers)  for  a  Single  Synodic  Cycle. 
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Solar  sail  trajectory  with  sail  profile 


Figure  35  Single  Cycle  Path  of  Solar  Sail  Cycler  with  Minimum  VJmars. 

A  noticeable  difference  between  the  solar  sail  cycler  and  the  traditional 
impulse  rocket  Aldrin  cycler  is  the  large  swingby  angular  deflections  with  respect  to 
Mars  and  Earth.  The  Aldrin  cycler,  because  it  is  minimizing  fuel  mass,  resides  in  a 
natural  Keplarian  orbit  most  of  the  time.  As  such,  it  tends  to  have  a  large  in  excess  of 
6  km/s  at  Mars  and  excess  of  5  km/s  at  Earth,  thus  restricting  turn  angles.  The  solar  sail, 
on  the  other  hand,  can  change  orbital  energy  with  no  impact  to  the  cost  function  and 
achieve  low  hyperbolic  excess  speeds  that  permit  large  turn  angles.  The  results  of  this 
analysis  show  that  a  75°  Mars  turn  angle  and  a  29°  Earth  turn  angle  provide  the  optimum 
path.  Interestingly,  the  sail  never  goes  down  to  the  minimum  allowable  pass  distance 
with  Mars  to  get  a  bigger  swingby  deflection.  Furthermore,  the  sail  swings  by  Earth  at 
the  maximum  allowable  perigee  distance,  not  the  minimum  allowable  distance.  Table  5 
summarizes  the  cycler  parameter  data  for  the  cost  function  with  y=  1 . 
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W  (kin/s) 

rp  ( planet  radii ) 

8 

Time  to  planet 

Mars 

2.53 

1.27 

75° 

1.1  months 

Earth 

4.3 

10 

29° 

18  months 

Table  5  Earth-Mars  Solar  Sail  Cycler  Flyby  Data  ((3=.I7,  y=l) 


It  should  be  noted  that  the  parameters  shown  in  Table  5  are  highly 
sensitive  to  the  states  right  before  the  Mars  and  Earth  encounters.  These  numbers 
represent  the  optimal  parameters  given  the  approximated  state  and  control  history.  By 
slightly  modifying  the  approximations  using  more  node  points,  states  preceding  a  knot 
could  change  enough  to  generate  a  slightly  different  optimal  parameter  set. 

To  verify  the  solutions,  states  were  propagated  using  initial  conditions, 
gravity  assist  conditions  and  the  DIDO -generated  control  history  using  a  Matlab®  ODE 
solver  as  described  in  the  Validating  Solutions  section.  Propagated  states  are  shown 
passing  through  the  DIDO  output  markers  in  Figure  34  producing  the  path  in  Figure  36. 

au  2 


Figure  36  Propagated  Path  (line)  with  DIDO  State  Output  (dots)  for  Single 
Cycle.  Solution  uses  115  total  nodes  (45  before  and  70  after  the  interior  knot). 
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To  capture  the  effectiveness  of  the  cyclic  end  condition  constraint,  the 
initial  and  final  propagated  states  were  compared  for  a  single  cycle.  A  comparison 
summary  appears  in  Table  6  where  the  error  between  final  and  initial  conditions  of  the 
propagated  path  are  shown  with  states  given  as  radial  distance,  r ,  Mars  lead  angle,  £  , 
velocity  magnitude,  V  ,  and  zenith  angle,  \|/  (complement  of  flight  path  angle). 


Final  state 

Value 

Initial  state 

Value 

%  Error 

rf 

1.0053 

ro 

1.000 

0.53 

Cr 

.6525 

Co 

.6512 

0.20 

yt 

1.1249 

n 

1.1256 

0.06 

V/ 

92.38° 

Vo 

93.86° 

1.60 

Table  6  Cyclic  End  Condition  Errorsfor  Synodic  Cycler 


Up  to  this  point,  we  have  only  minimized  the  hyperbolic  excess  speed  at 
Mars  with  the  weighting  factor  y  in  equation  (22)  equal  to  unity.  To  suit  the  needs  of 
any  particular  cycler  mission,  however,  it  may  be  desirable  to  minimize  a  combination  of 
at  both  Earth  and  at  Mars.  Varying  the  weighting  factor  produces  the  range  of  V^s 
shown  in  the  graph  in  Figure  37.  Ay  of  approximately  0.3  will  minimize  the  cost 

function  the  most  with  least  total  V». 


V-infinity  at  Earth  and  Mars 


'V-inf  Mars  V-inf  Earth 


Figure  37  Varying  y  in  the  Complex  Combination  Cost  Function 
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d. 


Remarks 


Using  a  reasonably  high  performance  solar  sail  to  achieve  optimal 
heliocentric  cycler  trajectories  is  a  viable  alternative  to  traditional  impulse  rocket  cyclers. 
There  are  several  interesting  ways  to  pose  an  “optimal”  solar  sail  cycler  problem.  One 
such  problem  would  be  a  synodic  cycler  that  achieves  Mars  and  Earth  encounters  with 
the  constraint  that  V,=0.  This  double  rendezvous  synodic  cycler  could  also  pose  an 
intriguing  design  optimization  problem  in  which  one  desires  the  minimum  sail  lightness 
number  required  to  achieve  a  double  rendezvous  in  a  synodic  period.  These  cases  are 
investigated  in  the  next  section. 

4.  Fun  with  Cycler  Trajectories 

a.  Double  Rendezvous  Synodic  Cycler 

Using  the  standard  performance  sail  ( (3  = .  1 7 ),  a  cycler  has  been 
presented  that  minimizes  V„  during  planetary  encounters.  It  may  be  interesting  to  seek  a 
synodic  cycler  that  is  constrained  to  rendezvous  with  each  planet  such  that  Vryj  is  zero  at 
both  planets.  Recall  that  performing  an  EME  double  rendezvous  mission  with  the 
standard  sail  with  unrestricted  end  conditions  took  at  least  2.41  years.  This  sail  would 
never  reach  Earth  again  within  a  synodic  period  of  2.135  years.  A  higher  performance 
sail  however  might  be  able  to.  Intuition  tells  us  that  there  ought  to  be  a  sail  with  just 
barely  enough  area  to  achieve  an  EME  double  rendezvous  within  a  synodic  period.  This 
sail  design  optimization  problem  may  be  stated  as  follows. 

Find  the  minimum  sail  lightness  number  that  would  enable  a  double 
rendezvous  cycler  where  V,*,  =  0  at  both  Earth  and  Mars. 

Minimize:  J  =  (3 

Subject  to  dynamic  constraints 

x-f(x,u)=0 
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and  event  constraints  given  in  Table  4  with  the  additional  constraint  that  the  end 
conditions  are  cyclic.  Hyperbolic  excess  velocity  is  zero  at  the  initial  and  final  times,  so 
the  only  additional  constraint  is  that  the  EM  lead  angle,  £  ,  is  the  same  at  both  times. 

C„=C, 

The  solar  sail  lightness  number  is  established  as  a  static  parameter,  a  sail 
design  characteristic  that  remains  constant  in  time.  The  results  in  Figure  38  and  Figure 
39  show  that  the  path  has  characteristics  of  both  a  synodic  cycler  and  a  double 
rendezvous.  There  is  symmetry  in  the  r  and  v  states,  and  antisymmetry  in  the  u  state  and 
control  angle  like  the  double  rendezvous.  Initial  and  final  conditions  are  precisely  the 
same  as  in  a  cycler,  although  no  gravity  assists  are  used.  The  time  of  flight  turns  out  to 
be,  of  course,  the  EM  synodic  period.  The  minimum  sail  lightness  number  required  to 
perform  such  a  mission  is  P=.297.  This  corresponds  to  a  sail  that  is  nearly  double  the 
size  of  our  standard  solar  sail  with  the  same  payload  mass! 


Minimum  Sail  Lightness  Number  for  EME  Double  Rendezvous 


Time  [years] 

Figure  38  Minimum  [3  Solar  Sail  States  and  Controls  for  an  EME  Double 

Rendezvous 
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Solar  sail  trajectory  with  sail  profile  history. 
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Lightness  no  =  .3  (min  beta  to  achieve  EME  double  rendezvous) 


Figure  39  Trajectory  for  Minimum  [3  EME  Double  Rendezvous 

b.  Taxi  Propellant  Cost 

One  mission  that  could  benefit  from  a  cycling  orbit  is  replenishing 
supplies  at  a  station  on  or  around  Mars.  In  concept,  a  taxi  craft  could  leave  its  parking 
orbit  about  Mars  and  greet  the  passing  cycler  sail  on  its  hyperbolic  trajectory  around 
Mars.  A  better  cost  function  in  this  case  would  be  the  fuel  required  to  meet  the  cycler  on 
its  swingby  path  from  a  parking  orbit  (equivalently,  we  could  minimize  the  A;’)-  We 
assume  a  circular  parking  orbit  with  a  radius  equal  to  the  closest  cycler  approach 
distance.  Since  the  same  exact  pass  conditions  would  be  met  every  cycle,  it  is  presumed 
that  the  taxi  craft  is  positioned  in  this  orbit.  F  ollowing  the  discussion  in  ref  16  p.  102,  the 
Av  required  to  go  from  a  circular  orbit  to  match  the  cycler’s  hyperbolic  orbit  is 


where  r  and  v  are  the  periapse  position  and  velocity  with  respect  to  Mars.  Since  the 

pm  pm  1  1  1  *'  1 

energy  of  a  hyperbolic  orbit  is 
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® hyp 


peri 

2 


and  the  periapse  velocity  of  the  spacecraft  with  respect  to  Mars  is 


v 


pm 


r 

pm 


The  change  in  velocity  required  to  join  the  cycler  craft  on  its  hyperbolic 
trajectory  is  then 


Note  that  if  the  circular  parking  orbit  grazes  the  closest  approach  point  of 
the  cycler  orbit  (for  fixed  r  )  the  Av  is  proportional  to  .  Thus  for  this  taxi  model,  a 

minimum  vmm  is  equivalent  to  a  minimum  taxi  Av  problem.  For  elliptic  or  3D  target 
orbit  models  where  the  closest  cycler  pass  distances  vary  from  visit  to  visit,  the  optimal 
taxi  intercept  path  would  differ.  In  this  case,  the  cost  would  also  have  to  include 
propellant  expended  to  get  from  a  nominal  parking  orbit  to  a  non -grazing  hyperbolic 
cycler  path. 


c.  Profiles  Using  Different  Sail  Performances 

As  with  any  other  blossoming  technology,  advances  in  sail  material  design 
are  related  to  the  amount  of  interest  and  thus  funding.  In  order  to  make  a  cycler  mission 
feasible  it  is  useful  to  know  what  mission  designs  are  available  for  any  given  sail 
performance.  In  the  next  analysis,  a  range  of  sail  lightness  numbers  were  fed  to  the 
dynamics  model  in  equation  (2)  and  analyzed  in  the  same  fashion  pres  ented  using  the 
standard  sail.  Results  for  y  =  1  in  the  cost  function  of  equation  (22)  for  a  range  of 
lightness  numbers  are  shown  in  Figure  40.  The  time  to  Mars,  t  ,  and  V  I  are  used  to 
compare  trajectory  characteristics.  The  higher  performance  sails  tend  to  take  a  longer 
time  to  reach  Mars  in  order  to  reduce  the  hyperbolic  exc  ess  speed  at  Mars. 
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Optimal  Cyclers  with  Different  Sail  Lightness  Numbers 
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Figure  40  The  Effect  of  Varying  Sail  Performances  on  a  Cycler  with  y  =  1  in  the 

Cost  Function 

As  the  sail  performance  is  reduced  down  to  (3  =  .01,  the  cycler  flyby 
characteristics  look  more  like  a  ballistic  Aldrin  cycler  that  passes  closer  to  Earth  using 
more  bend  angle.  The  sail  attitude  only  changes  the  orbit  enough  to  prevent  the  sail  from 
dipping  below  the  minimum  pass  distance  restriction.  A  sail  with  1/10 th  the  performance 
of  the  standard  sail  (i.e.,  |3=.0 17)  will  reach  Mars  in  a  short  time,  get  a  small  bend,  then 
sweep  outside  Mars  orbital  path  to  return  to  Earth  close  enough  to  get  a  large  bend.  The 
planet  encounter  parameters  are  shown  in  Table  7. 


VT  (km/s) 

rp  (planet  radii ) 

5 

Time  to  planet 

Mars 

8.4 

2.9 

7° 

5.2  months 

Earth 

5.3 

2.2 

61° 

20.4  months 

Table  7  Solar  Sail  Cycler  Characteristics  for  Cycler  (|3=.017,  y=l). 
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E  VALIDATION  OF  SOLUTIONS 

The  primary  method  of  validation  was  accomplished  through  numerical 
propagation.  Using  the  same  dynamic  equations  and  initial  conditions  employed  by 
DIDO  along  with  the  control  history,  a  propagator  will  generate  a  trajectory  that  may  be 
compared  with  the  output  state  history.  A  match  of  resulting  states  indicates  that  the 
solar  sail  with  the  given  control  history  follows  a  feasible  path  conforming  to  physical 
laws.  The  control  profile  may  be  obtained  from  the  direct  DIDO  output,  or  derived  from 
the  costate  history  (for  OCPs  without  an  interior  knot)  and  the  sail  control  law.  The  latter 
generally  provides  “smoother”  control  for  interior  controls  (when  control  angle  is  not  at  a 
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limit,  i.e.  -—<  a  <  —).  A  comparison  of  DIDO  and  costate-derived  controls  were 
shown  in  Figure  8for  the  flyby  and  Figure  15  for  the  rendezvous  problem. 

Numerical  propagation  was  accomplished  using  the  Matlab  ®  ODE45  and  ODE23 
solvers.  Controls  had  to  be  interpolated  between  DIDO  node  points  in  order  to  produce 
an  accurate  sail  attitude  at  time  steps  generated  by  the  ODE  solver.  Generally,  a  cubic 
interpolation  served  well  while  a  spline  method  proved  inaccurate  in  regions  with 
concentrated  node  points,  i.e.  near  knots. 

Validations  of  the  benchmark  problems  are  shown  in  Figure  9  and  Figure  17. 
DIDO  paths  and  the  propagated  paths  closely  match  as  shown  in  Table  8  listing  the  mean 
squared  error  for  the  different  runs  (DIDO  dynamic  constraints  were  met  with  less  than 
10 5  for  both  runs). 


State 

Flyby 

Rendezvous 

r 

1.2393e-009 

1.0944e-006 

e 

6.3148e-009 

3.3749e-004 

U 

1.0952e-009 

3.4343e-006 

V 

1.0258e-009 

3.4093e-006 

Table  8  Mean  Square  Error  of  DIDO  States  Compared  with  Propagated 
States  Using  Matlab®  ODE45 


66 


Propagation  of  swingby  trajectories,  such  as  those  used  in  the  cycler  problems, 
became  more  of  an  art  than  a  science.  The  propagation  was  enacted  for  one  cycle  only 
with  an  intermediate  event  condition  defined  at  Mars’  orbital  radius.  To  simulate  the 
gravity  assist,  a  switch  was  used  in  the  propagator  to  add  a  AV  to  the  sail’s  velocity  prior 
to  the  encounter  equal  to  the  AV  determined  by  DIDO.  Because  of  the  shaip 
discontinuity  in  the  path  velocity  due  to  the  gravity  assist  at  the  knot,  the  subsequent 
trajectory  was  driven  by  where  the  assist  occurred.  Small  differences  in  where  the 
outbound  path  encountered  Mars  caused  the  DIDO  states  and  propagated  results  to  differ 
on  the  inbound  leg  Figure  4 1 .  The  discrepancy  lay  in  how  the  sail  was  being  controlled 
just  prior  to  the  interior  knot.  Output  controls  often  appear  “shaky”  near  a  knot.  Wit  hin 
the  propagator,  these  somewhat  erratic  controls  must  be  approximated  at  the  time  steps. 
If  the  time  steps  of  the  propagator  were  large  in  comparison  to  the  DIDO  controls  at  the 
LGL  points,  approximation  errors  resulted.  These  small  errors  are  enough  to  slightly 
change  the  Mars  swingby  event  location  causing  the  remainder  of  the  trajectory  to 
deviate  from  the  DIDO  solution.  The  discrepancy  is  resolved  by  forcing  the  propagator 
step  sizes  to  be  approximately  the  same  size  as  the  distance  between  DIDO  -generated 
LGL  points  near  the  knot  (Figure  42).  In  this  way,  better  approximations  are  made  near 
the  defined  points  producing  closer  matches  between  DIDO  and  propagated  paths.  This 
can  be  accomplished  by  either  imposing  a  maximum  step  size  on  the  propagator  or  by 
adjusting  the  number  of  nodes  used  by  DIDO  to  manage  node  spacing.  Better  control 
approximations  in  the  neighborhood  of  the  planetary  encounters  yielded  a  closer  match 
between  the  optimal  solution  states  and  the  propagated  states.  This,  however,  only 
confirms  that  the  output  really  does  match  a  propagated  solution  for  the  given  erratic 
control  behavior  just  before  the  hard  knot.  We  desire  a  solution  that  does  not  require  the 
sail  to  perform  radical  attitude  changes  near  the  Mars  encounter. 
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Figure  41  Propagated  Path  with  “Missed”  Mars  Swingby  Phasing.  Propagated 
states  use  controls  interpolated  at  time  steps  different  than  DIDOs  LGL  distributed 
time  steps  causing  differences  at  the  interior  hard  knot. 


Figure  42  DIDO  Control  Output  at  LGL  Node  Points  and  Interpolated  Controls 
used  in  the  ODE45  Propagation  Near  a  Knot.  Step  sizes  match  fairly  well. 
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The  controls  may  be  made  “smoother”  by  adding  “inertia”  to  the  controls;  i.e.  by 
limiting  the  rates.  The  control  then  becomes  the  time  rate  of  change  of  the  cone  angle. 

x  =  [r  0  ,u  ,v  pc  f 

where  the  control  is  ua  =d(t). 

Using  this  control  and  state  variable  convention  provided  inherently  stable  cone 
angle  control,  a  more  desirable  design  for  a  solar  sail  attitude  control  system. 
Additionally,  the  pitch  rate,  a ,  was  be  restricted  to  provide  the  “inertia”  to  the  sail 
without  having  to  switch  from  a  3  DOF  problem  to  a  4  DOF  problem.  Numerical 
solution  time  was  expected  to  shorten  as  well  since  the  hodograph  of  “control”  a  is 
convex  (see  Appendix  E). 
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IV.  DEVELOPING  THE  OPTIMAL  THREE-DIMENSIONAL 

CYCLER 


We  now  extend  the  boundary  conditions  and  dynamics  of  the  2-D  solar  sail 
trajectory  models  of  the  previous  chapter  to  include  a  third  dimension  and  use  the  same 
methodical  approach  to  solve  the  optimal  cycler.  The  definition  of  a  “cycle”  for  the  3 
dimensional  problem  however  becomes  more  complex  than  the  2-D  case.  Initial  and 
final  cycle  end  conditions  do  not  repeat  exactly  in  a  synodic  period  as  they  did  in  the  2-D 
circular  coplanar  model.  Planetary  orbit  inclination  and  eccentricity  make  conditio  ns 
such  that  Earth  and  Mars  do  not  repeat  their  relative  positions  with  each  other  but  about 
every  15  years. 

The  motivation  for  obtaining  a  3-D  optimal  control  profile  for  a  solar  sail  cycler  is 
not  only  to  increase  the  fidelity  of  the  model  but  also  to  compare  results  to  the  2-D  case 
to  learn  more  about  the  nature  of  optimal  cycler  orbits.  The  inclination  of  Mars’  orbit 
with  respect  to  the  ecliptic  plane  is  roughly  1.85°,  not  a  great  difference  from  the 
coplanar  case.  Elliptic  orbits  for  Earth  and  Mars  are  modeled,  so  optimal  trajectories  are 
similar  to  the  elliptic  coplanar  solutions.  We  perform  the  same  series  of  benchmark  3D 
solutions  to  rate  them  against  their  2D  circular  and  elliptic  counterparts. 

A.  THREE-DIMENSIONAL  MODELS 

1.  Sail  Model 

Our  sail  in  the  new  model  now  has  an  extra  degree  of  freedom  due  to  the  addition 
of  the  third  dimension.  With  most  conventional  engines,  the  controls  generally  have 
three  degrees  of  freedom,  one  for  each  of  three  dimensions.  However,  with  the  sail, 
thrust  magnitude  is  dependent  on  the  cone  angle  providing  a  constraint,  therefore  only 
two  degrees  of  freedom  are  required.  The  cone  angle  serves  well  as  one  of  the  control 
variables  since  the  thrust  magnitude  is  related  to  it.  Another  angle  is  required  to 
determine  where  on  the  cone  the  sail  normal  vector  lies.  It  is  convenient  to  define  a  clock 
angle  as  the  angle  between  the  projection  of  the  sail  normal  onto  a  plane  normal  to  the 
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sun-line  and  a  reference  direction  in  the  same  plane;  so  we  adopt  a  model  used  in  ref  9 
p.  1 15.  This  reference  direction  is  taken  to  be  the  vector  normal  to  the  instantaneous 
orbital  plane  of  the  sail.  Figure  43  makes  the  representation  more  clear  with  unit  vector 
p  as  the  reference  direction.  Note  that  a  positive  clock  angle  rotation  is  in  the  negative  r 
direction  using  the  right  hand  rule. 


a  =  cone  angle 
8=  clock  angle 


V 


Figure  43  Solar  Sail  Control  Model  for  3  Dimensional  Dynamics 
2.  Dynamics  Model 

Since  the  thrust  magnitude  is  dependent  on  the  radial  distance  from  the  sun, 
spherical  equations  of  motion  provide  a  simple  way  of  including  the  acceleration  due  to 
the  sail  (see  Coordinate  Systems  Section) 

For  the  three  degree  of  freedom  state  variable  x  =  [r  B  ()  ,vr,(0e  ,co  |7  the  state 


derivative  is  given  by 
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where  the  last  three  terms  are  from  the  equations  of  motion  (Appendix  C). 

In  choosing  these  dynamic  equations,  we  must  bound  the  states  to  avoid 
singularities  at  r  =  0  and  cost])  =  0 .  It  turns  out  that  avoiding  these  state  values  keeps  the 
Jacobian  of  f  free  of  singularities  as  well. 


3.  Events  Model 

The  3D  target  orbit  manifolds  are  the  locations  and  corresponding  velocities  of 
Earth  and  Mars  along  their  respective  inclined  elliptic  orbits.  These  events  define  the 
boundary  conditions  of  the  optimal  control  problem.  In  defining  the  boundary 
conditions,  it  is  necessary  to  know  the  relative  orbital  shapes  and  orientations  of  the 
departure  and  destination  planets.  The  orbital  elements  of  Earth  and  Mars  orbits  are 
summarized  in  the  following  table 1 3. 


a  [AU] 

e 

i 

O 

to 

Earth 

1.0 

.0167 

0 

undefined 

102.9° 

Mars 

1.524 

.0935 

1.85° 

49.57° 

286.5° 

Table  9  Earth  and  Mars  Orbital  Parameters 
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In  all  the  problems  presented,  the  planet’s  true  anomaly  at  the  events  remains  a 
free  variable  to  be  determined  in  the  numerical  solution. 

In  order  to  rendezvous  with  a  planet,  it  is  essential  to  describe  the  position  and 
velocity  of  both  spacecraft  and  planet  in  a  common  3-D  (ND=3)  coordinate  system. 
Although  the  spacecraft  states  are  represented  in  spherical  coordinates,  the  events  at  t0 

and  tf  are  given  in  perifocal  coordinates  in  the  respective  planetary  orbital  planes.  To 

match  the  sail  states  with  the  planetary  states  at  the  end  conditions.  Earth  and  Mars 
orbital  planes  were  transformed  into  a  common  heliocentric  -ecliptic  coordinate  system 
along  with  their  respective  velocities  along  the  orbital  paths.  For  simplicity,  the  frame  is 
referred  to  as  the  E- frame.  These  orbital  states  define  the  target  manifolds  to  which  the 
spacecraft  event  conditions  are  to  be  constrained  (see  Appendix  D). 

First,  the  initial  Earth  angular  position  and  final  Mars  angular  position  are 
constrained  to  be  equal  to  the  sail’s  position  in  the  common  E-frame.  Resolving  these 
manifolds  in  the  Eframe  we  obtain  the  planetary  positions  in  a  common  coordinate 
system. 

Constraining  the  initial  and  final  sail  angular  positions  to  be  coincident  with  the 
initial  Earth  and  final  Mars  angular  positions  respectively,  we  obtain  for  small 
inclinations  (ref  4  p.  135)  in  the  E-frame 

9f0=eo-(®f  +^\)  9„,f  -0/-(com  +  n  ) 

where  0O  is  the  initial  spacecraft  position  and  0f  is  the  final  sail  position. 

In  a  perifocal  system,  the  planetary  polar  coordinates  rp  and  0P  are  related  by 

r  Pe  r  Pm 

l  +  eecos0eO  ""  l  +  e„,cos0m/ 

where  p  in  the  numerator  is  the  “parameter”  or  semi-latus  rectum  of  the  respective  orbit. 
Now  the  initial  and  final  sail  target  positions  in  Cartesian  coordinates  are  given  by 
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for  Mars 


're  o  cos(0o)^ 

'rni  cos(9/)> 

rosin(0o) 

l  0  , 

for  Earth  and  R  niE 

=  A 

m 

rnf  Sin(9/  ) 

0 

v  y 

where  Ae  and  Am  are  the  respective  Earth  and  Mars  3-1-3  rotations  transforming  into 
the  E- frame  (Appendix  D).  In  like  manner  the  velocities  of  the  planets  may  be  defined  in 
their  own  perifocal  frame  by 


where  V p  is  the  speed  of  the  planet,  |rp  |  is  the  distance  of  the  planet  from  the  sun,  and  ap 

is  the  semi-major  axis  of  the  planet.  Components  of  the  planetary  velocity  vector  are 
defined  in  terms  of  the  flight  path  angle  (1  by 


cos  (3  = 


r  \V 


and  sin  P/;  = 


h  V 

P  P 


'(e„) 


where  hp  represents  the  magnitude  of  the  planet’s  orbital  angular  momentum  vector. 

A  planet’s  velocity  vector  resolved  in  the  local  vertical  local  horizontal  frame 
(LVLH)  has  components 

u  =V  sin  P 

P  P  r  P 

vp  =V  cosp;, 

Finally  resolving  these  components  in  Cartesian  coordinates  and  transforming 
them  into  the  E- frame  yields  the  target  velocities  VeE  and  Vnin  (transformation  is  given  in 
Appendix  D).  For  the  rendezvous  mission,  V„f,  is  required  to  match  the  final  sail 
velocity,  but  knowledge  of  Mars’  velocity  is  not  needed  for  the  flyby  mission. 

The  constraints  may  now  be  set  such  that  the  spacecraft  states  at  the  initial  and 
final  events  are  equal  to  the  target  planetary  states.  The  initial  and  final  sail  positions  are 
given  as  follows. 
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x0  xf 


rse=  y0  yf 

Zq  Zf 

V  1  ) 

The  left  column  represents  the  sail  position  at  the  beginning,  r0 ,  and  the  right 
column  represents  the  sail  position  at  the  final  time,  r(  .  All  coordinates  are  resolved  in 
the  E- frame.  After  applying  the  transformation  given  in  equation  (45).  the  spacecraft 
velocities  at  events  in  LVLH  frame  are  given  by 


The  velocities  need  to  be  transformed  from  LVLH  into  Cartesian  coordinates  in 
the  E- frame  using  a  transformation  matrix  B  defined  as 

B  =  [Z?]=[tf3(-e)p2(f)] 

2— »3 


where  the  R  matrices  correspond  to  standard  rotation  matrices  (Appendix  D).  The 
spacecraft  velocity  in  the  Cartesian  E- frame  is  therefore 

VlB=  BV% 

The  left  column  of  \sE  represents  the  initial  sail  velocity,  V0,  and  the  right 
column  of  V  £  represents  the  final  sail  velocity,  V  .  Having  defined  all  the  necessary 

states  in  the  same  coordinate  sys  tem,  it  is  now  a  simple  matter  to  set  boundary  conditions 
for  the  initial  and  final  rendezvous  events. 


The  sail  position  starts  at  Earth  and  ends  at  Mars 


and  likewise  the  initial  sail  velocity  is  the  orbital  velocity  of  Earth  and  final  velocity  is 
equal  to  the  orbital  velocity  of  Mars. 
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(26) 


%-v£=o 


(27) 


v(-y„,=o 


The  final  velocity  event  constraint  in  equation  (27)  applies  to  the  rendezvous 
mission  only  where  the  spacecraft  must  match  the  velocity  of  Mars  in  its  orbit.  Table  10 
summarizes  the  event  conditions  in  equations  (24)  through  (27)  for  a  rendezvous.  All 
coordinates  are  resolved  in  Cartesian  coordinates  in  the  E-frame.  Within  the  DIDO 
framework  we  code  the  event  conditions  in  a  more  compact  form  that  ensures  end 
condition  manifolds  are  constrained. 

R  .  -  R  =  0  and  V  -V  =0 

st  et  st  et 

For  only  two  end  conditions  without  intermediate  conditions  the  dimensions  of  all 
the  matrices  are  ND  xNE,  where  ND  is  the  number  of  dimensions  and  NE  is  the  number 
of  events. 


Event 

Flyby  constraints 
[lower, upper] 

Rendezvous  constraints 
[lower, upper] 

Ro 

[R*>R*] 

Vo 

[V£,V£] 

[y£,vj 

R  / 

[R^b’R^] 

[RmE’R„£] 

V/ 

free 

[V£,v£] 

Table  10  Event  Conditions  for  3-D  Flyby  and  Rendezvous  Missions. 


B.  THE  OPTIMAL  CONTROL  PROBLEM 

In  keeping  with  the  step-by-step  approach  to  reaching  a  cycler  model,  we  first 
seek  the  minimum  time  solutions  to  the  3-D  flyby  and  rendezvous  missions.  To  this  end 
the  new  3D  OCP  is  stated  in  a  similar  fashion  as  the  2D  OCP.  Because  the  state  variable 
has  now  expanded  to  six  variables  (3  position,  3  velocity),  we  constrain  the  state 
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derivatives  to  the  3D  equations  of  motion  (equation  (23)).  States  are  bounded  to  avoid 
the  additional  singularity  at  the  solar  poles. 

C.  SOLAR  SAIL  CONTROL  LAW 

Since  the  sail  attitude  control  system  is  responsible  for  driving  both  the  cone 
angle,  a,  and  clock  angle,  8,  we  seek  two  steering  laws.  As  with  the  2-D  model,  we  turn 
to  the  principles  of  optimal  control  theory  to  generate  these  steering  laws.  Tailoring  the 
Bolza  cost  function  in  Lagrange  form  for  a  minimum  time  problem  we  get 

'/ 

J  =  E  +  J  Fdt 

h 

E  =  0,F=l 

Using  states  consistent  with  a  spherical  coordinate  system, 
x  =  [  r,0  ,<\>  ,vr,co0  ,ot  |7 ,  we  construct  our  state  derivative  vector  f  from  the  equations  of 

motion  (Appendix  C). 


(28)f(x,oc,5) 


2  M-  ,  PB 3 , 


rtOe~  cos~(|)  +  rxo^  --^r  +  ^cos  a 
r  r~ 


1 


rcos(|) 


2©e  ( sin(f>  —  vr  cos  <|> )  +  cos'asinasinS 


-2v  +-^-cos2a  sinacosS 


A 


—  coe  2  sin  <f>  cost]) 


Thus,  the  Hamiltonian  is 


/f(x,u,A.)  =  F  +  A/f(x  ,u) 
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H  =  1  +  Arr  x-AgG  +  A,0+Avr  r0 2  cos2(|)  +;t|)2--^-+-!^!icos  ,a 

\  r2  y2 


+  — - —  20  (n |)  sin(|)  -rcos(|))+-!-^-cos2asinasin8 
rcos(|)  |_  r2 


Ay,  ( (  Bn  ,  'l  • , 

H — —  — 2ft|>  H — —cos- (X sin CX coso  -0 -sin(])cos(]) 


The  Lagrangian  is  therefore 


Ux(  t),  u(r),  IQ)  =  H(x(  t ),  u(r),  A( f) )  +  p / (t  )g( x(  r),  u(Y) ) 


or  simply, 


L-H+  |i5S  +  u(/a  where  a  e  — — and  |Xg  =0  if  8  is  unconstrained. 


However  in  the  actual  implementation  of  control  bounds,  8  is  limited  by 
71  71 

8  e  -  —  to  avoid  non-distinct  controls  where  two  (a, 8)  control  pairs  could  produce 

the  same  resulting  thrust  vector.  It  will  be  apparent  later  why  this  restriction  is  important. 

First  deriving  the  control  law  for  the  clock  angle  we  calculate 

dL  L  P|t  2  .  c  K>  P|i  2  .  c 

—  =  — - - -cos-  a  sm  a  cos  8  — - — —cos-  a  sin  a  sin  8  +  (i5  =  0 

38  rcos(|)  r  r  r- 

Dividing  by  -^-cos2a  sin  a ,  we  obtain  the  tangent  steering  law  for  the  clock 
r 
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angle  8  interior  controls  for  a  ^0,±  — . 


tanS  =  - — — -  ,  ps  =  0 
cost)) 


The  control  angle  8  has  no  effect  on  the  resulting  thrust  vector  when  the  sail  is 

(  71  ^ 

exposing  its  full  area  to  the  sun  or  when  it  is  exposing  no  area  a=0,±  —  .  Applying 
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the  KKT  conditions  (Appendix  A)  to  obtain  the  circumstances  when  the  clock  angle  is  at 
the  “stops”  we  get 

5=-|,(i6<0 
8  =  +  ^H* 

where  the  sign  of  the  co vector  determines  which  way  the  clock  angle  is  oriented. 

Next,  we  desire  the  steering  law  for  the  cone  angle,  a .  Applying  the  same  KKT 
conditions  to  cone  angle,  a,  we  obtain  similar  results.  The  “stops”  of  the  cone  angle 
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control  occur  at  —  —  and  —  when  the  control  dual  variable  meets  the  following 
conditions. 


oc  =  --,Ha^° 
oc  =  +^-,na>0 


The  “interior”  cone  control  is  defined  where  the  dual  variable  is  zero. 


n  n  n 
<a  <—  ,u  =0 
2  2  “ 


Continuing  with  the  Minimum  Principle,  we  obtain 

dL  ,  Bu  /  2  .  Bu  .  _ ,  „  .  2  ,  \ 

—  =  -Kvr — —  3cos~asma  h - - -sm§  -2cosa  sin’  a  +  cos  a 

3a  r1  V  ’  r cos<|>  r2  V  ’ 


+  '<t>  ^  cosS  (-2cosa  sin2a  +cos3a )  +  (xa  =0 
r  r 


Dividing  by cos3  a  produces 


(30)  -3A,vr  tan  a  H - - — sin5  (-2tan2a  +  l)  +  — ^-cosS  (-2tan2a  + 1)  +  (Xa  =0 

r  cost])  v  r  v  ' 
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where  a  5*  ±  — .  This  corresponds  to  the  interior  controls  only  when  the  cone  angle 

multiplier  is  zero,  |Xa  =  0 .  Using  equation  (29)  for  tanS,  we  form  an  expression  for  cos 8 
and  sin  8  as  a  function  of  the  states  and  costates. 


cos  8  =  ± 


K cos(^ 


4  cos^t) 


sin8  =  ± 


+^,e  COS  <)> )  + 


The  sign  ambiguity  is  resolved  when  we  take  into  account  the  control  bounds,  i.e. 


8  e 


n  n 


,  therefore  sin  5  >  0 .  Recall  that  equation  (29)  is  not  valid  at  a  =  0 . 
Rearranging  equation  (30)  into  standard  quadratic  polynomial  form  for  tana 
atan2a+i?tana+c  =  0,  ae 
where  the  coefficients  are  given  by  the  following: 


,—  ]  and  a  5*0 
2  2 


a  =  -2—^2 — sin8  —2— — cosS 
r  cost])  r 

b  =  -3X 


^'e  -sin8+^-cos8 


rcos(|) 


Solving  the  quadratic  yields  the  tangent  steering  law  for  the  cone  angle  interior 
controls. 


(31) 


tana  : 


—b±>J  b1  —  4a  c 


2  a 


for  ae 


7t  31  ^ 
~2' 2  ) 


and  a  5*0 


Following  numerical  solutions  to  simple  optimal  control  problems  (i.e.  no  interior 
knots),  DIDO  will  output  the  time  histories  of  costate  as  well  as  dual  variables  associated 
with  the  control  angles.  Steering  laws  in  equations  (29)  and  (31)  use  these  outputs  to 


81 


generate  control  profiles  that  may  be  compared  to  the  DIDO  control  history  output  for 
verification. 


D.  RESULTS 

1.  Benchmark  Problem  Solutions  with  Elliptic  Inclined  Orbits 

The  question  in  reviewing  the  same  series  of  benchmark  problems  with  elliptic 
inclined  orbits  is  the  following.  Does  adding  a  degree  of  freedom  in  the  dynamics  to 
meet  the  3D  event  conditions  enable  us  to  further  reduce  the  cost  from  the  elliptic 
coplanar  case  or  does  it  hurt  us?  Apparent  from  the  first  several  solutions  is  that  the  3 
DOF  model  minimizes  the  transit  time  to  a  smaller  value  than  the  circular  coplanar 
model,  but  produces  a  larger  trajectory  time  than  the  elliptic  coplanar  case. 

a.  Earth-Mars  Flyby 

To  reach  Mars  as  quickly  as  possible  and  sail  past  without  gravitational 
interaction,  it  was  intuitive  that  Mars  periapsis  provided  the  minimum  path  distance  in  the 
coplanar  model.  Now  when  we  consider  that  the  target  orbital  planes  are  inclined  with 
respect  to  each  other  it  may  be  desirable  to  choose  a  path  that  avoids  orbit  “cranking” 
even  when  hitting  the  periapsis  is  not  possible  in  plane.  The  solution  shows  that  the  sail 
changes  planes  a  small  amount  in  order  to  reach  Mars  at  its  perihelion  as  shown  in  Figure 
44.  Mars’  orbit  has  a  small  inclination  with  respect  to  the  ecliptic  so  it  is  preferable  to 
traverse  a  short  distance  even  though  a  small  cranking  maneuver  is  incurred.  Figure  45 
reveals  that  the  time  to  intercept  Mars  with  our  standard  sail  ( (3=.  17)  is  138  days,  which  is 
about  the  same  as  transiting  to  an  in -plane  elliptic  Mars  orbit  (137  days).  The  sail  cone 
angle  control  history  is  very  similar  to  the  coplanar  case  cone  history  while  the  clock 
angle  steadily  increases  (Figure  46).  Controls  derived  from  the  tangent  steering  law 
using  costates  are  compared  to  the  DIDO  control  output  showing  a  close  match. 
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Mars  flyby  sail  trajectory.  DIDO  states  and  propagation 


Figure  44  Flyby  Mission  to  Mars  with  Elliptic,  Inclined  Planetary  Orbital 
Planes.  DIDO  output  (dots)  and  propagated  path  (line).  Mars  orbit  inclination  is 
exaggerated  for  display  purposes. 


83 


Mars  flyby.  DIDO  states  and  propagated  states.  20  nodes. 


Figure  45  State  History  for  Mars  Flyby  (3D  Model) 


Steering  law  with  Costates  vs.  DIDO  output. 


Time  (years) 

Figure  46  Cone  and  Clock  Angle  Controls  for  Mars  Flyby.  History  is  shown  for 
DIDO  (markers)  and  tangent  steering  law  using  states  and  costates  (lines). 
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b.  Rendezvous 

Observing  the  (])  and  C0o  states  in  Figure  47  and  Figure  48  gives  insight 

into  what  is  required  to  perform  the  plane  change  maneuver.  The  sail  trajectory  remains 
in  the  Earth  ecliptic  plane  for  the  first  quarter  of  the  mission  time,  then  “cranks”  the 
orbital  plane  to  match  the  final  sail  (])  and  C0o  with  Mars  exactly  at  the  time  of  Mars 

encounter.  The  mission  is  1.01  years,  which  is  a  month  faster  that  the  circular  coplanar 
and  takes  only  marginally  longer  than  the  elliptic  coplanar  orbit  model  (0.977  years). 


80  nodes 


Time  (years) 

Figure  47  EM  Rendezvous  with  Elliptic,  Inclined  Orbits. 
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80  nodes 


Figure  48  Optimal  Profile  for  <j)  and  co^  States  in  the  3D  Rendezvous  Mission. 

Turning  to  the  costate  analysis,  it  is  evident  that  as  in  the  2-D  solutions 
there  is  a  critical  time  during  the  maneuver  at  which  all  costates  are  at  a  local  maximum 
or  minimum  or  zero  (Figure  49).  It  is  beyond  the  scope  of  this  thesis  as  to  why  this 
behavior  occurs,  however  at  this  point  in  time  (about  1/4 *  into  the  mission  time)  the 
control  angles  make  large  adjustments  to  achieve  the  state  changes  previously  described. 
When  inserted  into  the  control  laws  derived  in  equations  (29)  and  (31),  the  costates  and 
states  will  produce  derived  cone  and  clock  angle  controls.  A  comparison  of  costate- 
derived  controls  and  DIDO  controls  is  shown  in  Figure  50.  The  optimal  control  does  not 
rotate  the  cone  angle  to  the  extremes  that  the  coplanar  model  control  did.  The  controls 
here  only  pressed  the  cone  angle  between  18  and  60  degrees  versus  between  10  and  70 
degrees.  The  optimal  cone  angle  has  to  account  for  the  clock  angle  as  well  in  this  3  DOF 
model. 
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Costates  for  the  minimum  time  3-D  Mars  rendezvous 


Figure  49  Costate  History  for  Minimum  Time  Rendezvous  with  3D  Orbits 

Model. 

Checking  feasibility,  the  propagated  controls  were  observed  to  match 
closely  with  the  DIDO  output  (Figure  51).  For  optimality,  the  Hamiltonian  is  constant  in 
Figure  52. 
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Derived  controls  using  costates  and  DIDO  controls 


Figure  50  Control  History  for  3D  Rendezvous.  Displayed  are  DIDO  controls 
(markers)  and  the  controls  derived  from  the  tangent  steering  law  with  the  states  and 

costates  as  inputs  (lines). 


Propagated  path  and  DIDO  state  output 


Figure  51  Propagated  Path  of  3D  Rendezvous  ( line)  with  DIDO  Output  (dots). 

Mars  inclination  is  exaggerated  for  display  purposes. 
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Hamiltonian  for  3D  EM  rendezvous.  80  nodes 


Figure  52  Hamiltonian  of  3D  EM  Rendezvous  Problem 
c.  Double  Rendezvous 

To  optimize  the  Earth -Mars-Earth  (EME)  double  rendezvous,  we  add  to 
the  rendezvous  event  conditions  a  return  trip  to  Earth.  The  2-D  elliptic  coplanar  double - 
rendezvous  exhibited  symmetry  in  states  and  controls  on  the  outbound  and  inbound  legs 
(Figure  31).  Now  using  inclined  ellipses  with  their  corresponding  positions  and  velocity 
vector  fields  as  target  manifolds,  we  seek  insights  into  the  double  -rendezvous  with  this 
higher  fidelity  model. 

The  event  conditions  now  include  an  interior  event  corresponding  to  a 
rendezvous  with  Mars  and  a  final  event  corresponding  to  a  rendezvous  with  Earth.  To 
ensure  Earth  is  actually  present  when  the  sail  intersects  its  orbit  manifold,  we  constrain 
the  time  of  flight  using  Kepler’s  equation.  This  last  non-linear  constraint  replaces  the 
corresponding  one  in  the  2-D  model  that  was  listed  in  Table  4.  Fortunately,  the  non¬ 
linear  equation  does  not  need  to  be  solved  explicitly  as  it  simply  serves  as  a  constraint  to 
the  trajectory.  In  this  way,  the  time  of  flight  is  coupled  with  the  eccentric  anom  aly  that 
the  Earth  traverses,  E .  Earth’s  angular  displacement  traversed  is  constrained  by  the 
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flight  time  of  the  spacecraft.  As  long  as  Earth’s  true  anomaly  change  is  the  same 
displacement  as  the  spacecraft  angular  displacement  plus  a  multiple  of  271  phasing  will 
have  been  achieved. 


In  order  to  rendezvous  with  a  planet,  we  need  to  describe  the  position  and 
velocity  of  both  spacecraft  and  planet  in  3-D  (ND=3)  for  all  4  events  (NE=4). 

Spacecraft  3D  position  at  events  in  Eframe  Cartesian  coordinates  are 

expressed  as 


( 

*0 

■C 

xl  xf 

.Vo 

fi 

yt 

V/ 

<  z/  J 

velocity  at  events  in 

LVLH. 

f 

Co 

Vri 

Vri 

vs  = 

Veo 

vii 

Vi 

ve  / 

> 

V 

V.-  v 

Transforming  the  sail  velocity  vector  from  LVLH  into  in  E-frame,  we  use 
the  transformation  matrix. 

B  =  [B]=[tf3(-e)p2(4»] 

2— >3 

So  spacecraft  velocity  in  E- frame  Cartesian  coordinates  is 


V  =BV 

sE  s  sE 


Because  the  craft  is  returning  to  Earth,  it  must  be  properly  phased  with 
Earth  to  ensure  the  planet  is  there  when  the  spacecraft  arrives.  This  constraint  makes  an 
Earth  orbital  motion  model  necessary  accounting  for  its  position  in  space  and  time.  The 
following  equations  are  used  to  constrain  initial  and  final  Earth  events  to  ellipses  in  three 
dimensions.  All  state  vector  multiplication  and  division  operators  are  element-wise. 

First,  we  define  the  Earth  initial  and  final  angular  positions  in  the  perifocal 

frame. 
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—  [®eO  ]  —  [^0  ®/]  \,le  ,le  ]  » 

Notice  that  0e  e  R 2  where  [0O  0/l  are  the  initial  and  final  angular 

spacecraft  positions  in  E- frame  which  are  coincident  with  Earth’s  initial  and  final 
position.  The  following  equations  relate  Earth’s  eccentric  anomaly  to  its  true  anomaly. 
Eccentric  anomaly  is  expressed  as  a  vec  tor  with  each  element  corresponding  to  a  knot,  or 
end  condition.  The  purpose  of  expressing  it  in  this  manner  is  to  be  able  to  handle 
multiple  cycles  with  multiple  encounters  (  NE  >  1 )  using  a  single  variable. 


The  initial  and  final  cosine  of  Earth’s  eccentric  anomaly  are  contained  in 
the  vector  cosE.  The  vector  operations  are  performed  in  an  “element-wise”  fashion 
where  the  cosine  of  the  vector  0(  is  the  cosine  of  each  individual  element  of  vector  0, . 


e  +cos(0  ) 

cosE  =  -£ - ,  cosEe  RN* 

l  +  ee  cos  ( 0„ ) 

Earth’s  radial  position  as  a  function  of  cosE  is  given  by 
r  =a  [l-e  cosE)  ,  r  e  MW£ 

e  e  \  e  /  7  e 


where  r  =  [r  0  ,ref  ]  so  that  vector  sinE  may  be  expressed  as 


sinE  = 


r,  sin(  0e ) 


,  r  e 


ae^~ee 

Thus,  Earth’s  initial  and  final  eccentric  anomaly  are  contained  in  the 


vector  E . 


E  =  tan  1 


sinE 

cosE 


,  Eg: 


Earth’s  eccentric  anomaly 


The  next  set  of  equations  is  used  to  define  Earth’s  position  based  on  the 
initial  and  final  time. 


Earth’s  mean  motion 
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Initial  mean  anomaly 


M0=E0  +  eeSinE0 

M  =  netTRANSIT  Final  mean  anomaly 

Now  we  define  the  initial  and  final  Earth  position  event  matrix  in 
heliocentric  Cartesian  coordinates.  E-frame  where. 


with  velocity  components  expressed  as 

cos II,  =  — ,  cosB  6l"! 

e  r  ®  V 

|  e  |  e 

sinB  =  —  ersin(8  ),  sinB  el*1 

e  /  -tt  E  V  e  /  5  e 

k  V 

e  e 

where  cosBe  and  sinBe  are  the  cosine  and  sine  of  Earth’s  flight  path  angle  at  the  event 
knots.  Earth’s  velocity  in  radial  and  transverse  components  in  the  perifocal  frame  are 
expressed  in  the  following  equations  where  the  symbol  ®  is  used  to  emphasize  that  this 
is  an  “element-wise”  multiplication  such  that  each  element  in  one  vector  is  multiplied  by 
the  corresponding  element  in  the  other  vector. 

ue  =  Ve  ®sinl\  ,  ue  e  M"* 

v  =  V  ®cosB  ,  v  g  M"* 

ut 

VeCART=[^(-0e)]  ,  VeCARTGMW^ 

0 

V 
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VeE  e 

where  Ae  is  a  matrix  performing  a  3-1-3  transformation  from  Earth  perifocal  to  Cartesian 
heliocentric  inertial  (E- frame).  See  Appendix  D  for  the  derivation. 

Having  defined  all  the  necessary  quantities,  it  is  now  a  simple  matter  to  set 
boundary  conditions.  Letting  the  spacecraft  initial  and  final  event  states  be  RsEnd  and 
VsEnd,  we  have 


V  =  A  V 

T  eE  e’eCART 
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and  intermediate  event  conditions  are 
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so  the  boundary  conditions  at  the  earth  encounter  events  become 


and 


We  can  establish  the  Mars  intermediate  events  in  a  similar  fashion. 

R  -R  _ =0 

i  mt 

and 

v,-v„,£=  0 

One  last  event  constraint  ensures  proper  phasing  with  Earth.  This  is 
accomplished  by  forcing  Kepler’s  equation  as  a  constraint. 
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E-eK  sinE-AC  —M.  =0 

E  Of 

Earth  mean  anomaly  and  eccentric  anomaly  are  represented  by  M  and  E 
respectively.  The  event  constraints  for  the  3D  EME  rendezvous  mission  are  summarized 
in  the  table  below. 


Event 

Double  rendezvous 
constraints 
[lower, upper] 

End 

[vE,v,l 

R 

l 

t  tvE  ’  tnE  ^ 

v, 

[v^v,] 

E-er  sin  E  -  M  -  M  , 

E  Of 

[0,0] 

Table  11  EME  3D  Double  Rendezvous  Event  Constraints 

States  and  controls  for  the  3  DOF  model  are  bounded  by 

xeS  cl‘ 
ueU  czM2 

Where  the  set  S  is  chosen  to  avoid  singularities  in  the  dynamics  and  the  set 
U  was  chosen  to  avoid  duplicity  in  the  controls.  Bee  ause  it  was  desired  to  obtain  distinct 
(a  ,8  )  pairs  (every  sail  normal  vector  orientation  is  represented  by  only  one  (a  ,5  )  pair), 

7t 

the  bounds  were  such  that  0  <a  <  —  and  -n  <8  <7t  .  It  was  recognized  that  the  sin  and 

cos  components  of  the  normal  vector  could  have  been  used  as  the  controls  with  proper 
path  constraints,  however  the  numerical  solutions  took  longer  to  complete,  so  the  cone 
and  clock  angles  were  retained  as  the  controls. 

Results  and  Analysis 

Once  more,  Mars  aphelion  played  the  role  of  the  optimal  rendezvous  point 
at  the  intermediate  event  for  the  round  trip.  Modeling  Mars’  small  inclination  with  the 
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ecliptic  produced  a  solution  that  looked  very  similar  to  the  2D  ellip  tic  when  case  when 
looking  from  the  top  view  (compare  Figure  53  with  Figure  31).  The  mission  time, 
however,  was  slightly  longer  than  the  2D  case  taking  2.42  years  instead  of  2.3  years  since 
the  spacecraft  must  change  planes  twice.  The  mission  time  of  the  3D  model  is  the  same 
as  that  of  the  2D  circular  orbits  model  except  that  there  is  no  longer  the  precise  symmetry 
of  states  and  controls  exhibited  by  the  other  model.  In  three  dimensions,  the  spacecraft 
takes  1.03  years  on  the  outbound  leg  and  1.39  on  the  inbound  (2D  circular  coplanar 
model  required  1.205  years/leg).  One  of  the  amazing  features  of  the  path  is  that  the 
spacecraft  sails  out  of  both  Mars  and  Earth  orbital  planes  for  a  significant  portion  of  t  he 
return  trip  to  Earth  rendezvous  (Figure  54).  This  large  departure  from  the  planetary 
planes  is  accomplished  in  order  to  reach  Earth  orbit  in  phase  with  Earth  to  an  accuracy 
2.64E-4  radians  (-39E3  km  from  the  center  of  the  Earth).  Testing  the  feasibility,  the 
control  profile  was  entered  into  the  propagator  producing  the  result  in  Figure  55. 


Figure  53  3D  EME  Double  Rendezvous  Path  (top  view) 
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Figure  54  3D  EME  Double  Rendezvous  Path  (oblique  view) 


Figure  55  Propagated  Path  (line)  and  DIDO  States  (dots)  for  3D  EME  Double 

Rendezvous 
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2.  Earth-Mars  Cycler  Set-up 

Formulating  the  3-D  elliptic  orbit  cycler  model  into  an  OCP  is  a  formidable  task 
because  cyclic  end  conditions  are  no  longer  imposed  every  EME  round  trip.  Many  round 
trips  occur  before  Mars  and  Earth  positions  are  exactly  the  same  as  they  were  at  the  start 
of  the  trajectory.  The  events  as  structured  in  the  previous  section  are  set  up  to  handle 
multiple  passes  for  a  “cycle”.  Since  meeting  synodic  conditions  with  each  round  trip  is 
no  longer  required,  a  possible  OCP  formulation  would  be  to  seek  the  minimum  time  to 
achieve  a  given  number  of  EM  passes.  Using  the  framework  outlined  in  this  thesis,  this 
would  entail  performing  successive  solutions  to  increasing  number  of  EM  visits  (and  thus 
number  of  knots)  using  intercept  end  conditions  until  initial  state  conditions  are  achieved 
at  the  final  time. 
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V.  CONCLUSIONS  AND  FUTURE  WORK 


A.  CONCLUSIONS 

Solar  sails  provide  a  viable  means  of  propulsion  for  an  Earth-Mars  cycler.  For  a 
solar  sail  with  lightness  number  (3  =.17,  an  orbit  can  be  maintained  with  rapid  revisit 

times  between  destination  planets  and  slow  V^s  at  each  planet,  both  desirable  cycler 
characteristics.  Furthermore,  the  cycler  lifetime  is  only  dependent  on  sail  degradation  in 
the  space  environment  since  no  propellant  is  depleted  making  solar  sail  cyclers  attractive 
in  comparison  to  conventional  low  thrust  and  chemical  propellant  cyclers. 

The  framework  in  this  thesis  outlines  a  method  of  state  and  control  optimization 
for  both  2-D  and  3-D  cycler  models.  This  framework  can  be  applied  to  other  missions  as 
well  by  changing  the  event  conditions.  With  relatively  simple  dynamics  and  sail  models, 
trajectory  design  could  be  accomplished  by  properly  formulating  desired  events  in  the 
form  of  constraints.  The  pseudospectral  method  used  within  DIDO  to  solve  the  OCPs 
produces  optimal  controls  and  states  that  were  verified  as  feasible  (through  propagation 
of  initial  conditions  with  a  third  party  ODE  solver)  and  optimal  using  the  necessary 
conditions  for  optimality. 


B.  OPPORTUNITIES  FO  R  FUTURE  WORK 

The  following  are  some  areas  of  interest  for  continuing  work. 

Earth- Venus-Mars  cycler.  Is  a  cycler  possible  with  three  target  planets  and  if 
so,  what  does  it  look  like?  Is  it  possible  to  reduce  the  revisit  times  from  a  two -planet 
cycler?  This  problem  would  require  the  solution  to  include  the  optimal  order  of  planetary 
encounters. 

Use  a  solar  sail  to  establish  an  optimal  halo  orbit.  Many  non-Keplarian  orbits 
are  possible  with  solar  sails.  Orbits  around  Fagrange  points  may  be  modeled  and 
examined  for  missions  such  as  early  detection  of  geomagnetic  storms  on  the  sunward  side 
of  the  Earth-Sun  El  point  (ref  9  p.  231).  Cycling  orbits  between  L  points  could  also  be 
explored. 
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Implement  two  synodic  period  cyclers  with  solar  sails.  Reference  19  outlines  a 
plan  to  make  multiple  passes  of  a  target  planet  within  a  cycle  (2  synodic  periods).  Planet  - 
centered  dynamics  may  play  a  role  if  trajectories  spend  significant  time  in  a  planet’s  SOI. 

SEP/Solar  Sail  hybrid  cyclers.  Modeling  a  SEP/solar  sail  hybrid  spacecraft 
would  show  when  it  is  optimal  to  favor  the  SEP  and  when  it  is  preferable  to  favor  the  sail 
for  various  cycler  orbits  (or  any  trajectory  for  that  matter). 

Variable  Isp  cycler.  When  is  it  optimal  to  use  high  thrust-low  Isp  and  when  is  it 
best  to  use  low  thrust-high  Isp  in  a  cycler  trajectory?  There  are  trades  to  examine  with 
weighted  cost  functions  between  fuel  mass  and  time. 

Make  smoother  controls  near  hard  knots .  Changing  the  controls  to  include 
“inertia”  assisted  in  making  controls  somewhat  smoother  near  knots,  however  for  3D 
problems  it  was  important  to  implement  constraints  to  ensure  that  arbitrary  sail 
orientations  resulted  from  unique  cone  angle/clock  angle  pairs.  Implementing  t  he  sin  and 
cos  components  of  these  angles  as  controls  (with  a  path  constraint)  slowed  the  algorithm 
down  and  was  not  implemented,  however  this  may  have  alleviated  some  of  the  jumpy 
control  outputs. 

Improve  gravity  assist  model  The  matched  asymptotes  model  as  presented  in 
this  thesis  is  accurate  for  the  problems  posed,  however  more  accurate  solutions  to  some 
problems  may  call  for  a  better  gravity  assist  model.  Possibilities  are  to  either  change 
dynamics  inside  SOI  or  use  a  jump  in  time  and  theta  at  the  in  terior  knot  in  addition  to 
velocity. 

Improve  sail  model.  A  higher  fidelity  model  of  the  sail  is  available  that  includes 
billowing  and  diffuse  reflection  (see  ref  9  pp.51-54).  Once  the  dynamics  and  events 
models  have  been  thoroughly  tested,  there  are  some  design  trades  that  could  be  explored 
as  far  as  the  sail  itself. 

Change  to  dynamic  equations  that  are  better  suited  for  cyclers  .  The 

discussion  on  coordinate  systems  outlined  some  benefits  and  pitfalls  of  using  certain 
coordinates  in  cycler  trajectories.  Perhaps  other  coordinates  than  those  presented  here 
can  be  better  implemented  for  optimal  cyclers  (e.g.  equinoctial). 
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Solar  sail  stability.  Solar  sails  have  been  used  in  spacecraft  to  dump 
accumulated  momentum  in  reaction  wheels  due  to  torque  perturbations.  A  possible  OCP 
would  be  to  minimize  the  time  to  dump  momentum  in  reaction  wheels  for  a  given  earth  - 
centered  orbit. 
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APPENDIX  A.  APPLICATION  OF  THE  MINIMUM  PRINCIPLE 


10,17 


The  first  step  in  solving  an  optimal  control  problem  is  to  construct  the  scalar 
Hamiltonian  function,  H , 

H  (x,u,  t.  A.)  =  F  (x,u,t)  +  X(t)T { (x,u,  f) 

where  f(x,u,f)  are  the  dynamic  constraints  on  the  system,  and  X(i)  are  the  Lagrange 

multipliers  called  costates.  According  to  Pontryagin’s  Minimum  Principle,  the  optimal 
control  is  obtained  by  solving  the  following  problem  at  each  instant  in  time. 

Minimize  H  with  respect  to  u,  with  u  e  U  where  U  is  the  set  of  all  allowable 
controls.  To  solve  this  problem,  the  Lagrangian  of  the  Hamiltonian  is  written  as: 

L(x,u,f,A,,ns)  =  //(x,u,  f,A,)+ps  (f)Tg(x,u,f) 

where  g(x,u,t)  are  the  path  constraints  and  p  (t)  are  the  associated  path  covectors. 

The  path  constraints  include  all  trajectory  path  constraints  as  well  as  any  state  and  control 
bounds.  Applying  the  Karush-Kuhn-Tucker  (KKT)  theorem  to  the  Lagrangian  results  in 
a  set  of  first  order  necessary  conditions  and  provides  a  means  to  demonstrate  the 
optimality  of  a  solution. 


(32) 


P*  (f)Tg(x,u,r)=0 

with 

^0  g,  =  g(x,u,x) 

(  =  >0  g(x,u,x)=g„ 

=0  '  g,  <g(x,u,x)<g„ 

any  g,  =  gu 
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The  third  case  above  describes  a  special  condition  when  the  constraints  in  g  are 
interior  or  non-binding, 

g,  <g(x,u,f)<g„ 

For  these  cases,  the  multipliers,  p  =0  and  equation  (32)  simplifies  to 


c)L  _  c)H  _ 

3u  3u 

A  necessary  second  order  condition  for  optimality  when  p  =  0  is  that 


(33) 


d2H 

3u2 
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APPENDIX  B.  DERIVATION  OF  THE  2-D  EQUATIONS  OF 

MOTION 


Lagrangian  mechanics  are  used  to  achieve  the  2-D  equations  of  motion  using 
polar  coordinates.  Velocity  variables  are  then  written  into  a  more  familiar  form  using  the 
radial  and  transverse  velocity  components.  Figure  1  shows  the  states  and  coordinate 
convention  used  in  the  derivation.  In  general  form  Lagrange’s  equation  is  written  as 


d  (dL  8 
dt  ydx 


where  L  i  the  Lagrangian,  x  is  the  generalized  coordinate  vector,  and  Qnc  is  the 
generalized  non -conservative  force  vector18. 

We  start  with  writing  the  specific  kinetic  energy  as 


T  =—  v-  v  =  —  (r2  +  r202 ) 
2  2 


Potential  energy  due  to  gravity  is 

v_  GM _  p 
r  r 

The  Lagrangian  may  be  written  as 

L=T-V  =-(f2  +r2d2)  +— 
2  r 


The  non-conservative  generalized  force  is  due  to  the  solar  radiation  pressure 
(SRP)  on  the  solar  sail.  Acceleration  magnitude  is  given  by 


T_ 

m 


Thus  the  acceleration  in  the  radial  direction  is 


T  pp  3 

a,  = — cosa  = — j-cos  a 
m  r0~ 


and  the  applied  acceleration  in  the  along -track  direction  is 
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T  .  P(I  2  . 

ae  = — sina  = — -cos~  a.  sin  a 
m  r0 


Writing  the  Lagrangian  equations  of  motion,  we  obtain 


d 

dt 


rdL'' 

\drJ 


=  r 


—  =  r62--^- 

dr  r1 


(34) 


(35) 


r-rG2+  -^-=  a 


dt 


f  0L  ^ 

vae  y 


:  2rid  +  r2©  —  =  0 
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2rrQ  +r2Q 


■  aar 


where  aBr  is  a  generalized  specific  torque. 

Equation  (35)  reduces  to  the  conservation  of  angular  momentum  when  aB  —  0 . 
Now  to  convert  the  r  s  and  0  s  into  u  and  v  velocity  components,  we  use  the 
following  relations. 

Transverse  velocity  is 
(36)  v  =  r0 

and  transverse  acceleration  is 


(37)  v  =  r0  +  rO 

Using  the  relation  in  equation  (36)  with  equation  (34)  we  write  the  equation  of 
motion  for  u . 


(38) 


Pi*  3 

,  +  4T-  cos  a 
r  rn 


Rewriting  equation  (35)  and  factoring  an  r  yields 
(39)  r  (r0  +  rQ  +  r0 )  =  ae  r 
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Recall  that  u  =  r  therefore  uv  =  iiQ  .  Using  (37)  with  ( 39)  we  obtain  the  equation 
of  motion  for  v . 

rv  +  rrd  =  aBr 
uv 

v  = - +  ae 

r 

(40) 
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APPENDIX  C.  DERIVATION  OF  THE  3-D  EQUATIONS  OF 

MOTION 


As  with  the  2-D  case,  we  use  the  Lagrangian  approach  to  acquire  the  3-D  equations  of 

motion. 


Figure  56  shows  the  RSW  spherical  coordinate  system  (ref  4  p.  42)  and 
convention  used  to  represent  the  states. 


z 


Figure  56  RSW  Coordinate  System 
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In  general  form  Lagrange’s  equation  is  written  as 


d  '  dL  'j  dL  T 
~dt\dx  j~fa~Qnc 


where  L  is  the  Lagrangian,  x  is  the  generalized  coordinate  vector,  and  Qnc  is  the 
generalized  non-conservative  force  vector18.  The  generalized  non-conservative  force  for 
this  case  is  the  solar  radiation  pressure  (SRP)  on  the  solar  sail  resolved  in  radial  ( R ), 
along-track  ( S ),  and  cross-track  ( W )  directions. 

We  start  by  writing  the  specific  kinetic  energy  as 


where  the  velocity  v  is 


r=-vv 

2 


v  =  rR  +  ?e  cos  <]xS  +  r(|)W 


Specific  potential  energy  is 


V  =  -  — 


making  the  Lagrangian 


L=T -V = 


V  =  —  r 2  +  (r9  cos([) )  +  (r(f>  j 


The  external  acceleration  due  to  the  SRP  on  the  sail  has  a  magnitude  of 

T  pp  2 
a  =  — =-t-^-cos~a 


and  have  components  in  the  radial,  along-track  and  cross- track  directions,  written  as 
a  cos  a  R  +  a  sin  a  sin  8  S+  a  sin  a  cos  8  W 

Having  established  the  necessary  quantities,  we  proceed  to  derive  the  Lagrangian 
equations  of  motion.  For  the  r  state,  we  have  the  following. 


d  (dL 
dt  dr 


—  =  r6 2  cos2 (|)  +  n|) 2  - -^- 
dr  r~ 
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(42) 

Turning  to  the  0  state  we  form  the  quantities 

\  ... 

—  — r  =  r20  (-2<])  cost))  sin())  )  +  cos2t|)  (r2©  +  2r?0  )  —  =  0 

dt  ^80  j  ((0 

The  generalized  force  is  a  torque  about  the  W  axis 
Q] v  =t/sinasin8  ( r cos cf) ) 

Lagrange’s  equation  (41)  becomes 

r2©  ( — 34>  cost])  sint])  )  +  cos2t|)  (r2©  +  2r?0  j  =  a  sin  a  sin  8  (r  cost]) ) 

Dividing  by  r  cost])  results  in 

?0  ( — 2sin t|)t]) )+  cost]))0  +  2f0  =  asina  sin8 
Finally,  we  rearrange  terms  to  get  the  equation  of  motion  for  0  . 

(43) 

Note  that  we  need  to  be  cautious  at  the  singularities  where  r=0  and  cos  t])=0. 


Now  for  the  tf)  state  we  write 


d  f dL^ 
dt  d§  / 


=  2rrtf>  +  r2tf> 


—  =  -r0 2  cost])  sintb 
3t|) 


The  generalized  force  in  the  t|)  sense  direction  (torque  in  this  case)  is  given  by 

Q ^  =  arsina  cosS 

Substituting  into  equation  (41)  yields 

2nt\>  +r2t])  +  r202  cost])  sint])  =arsinotcosS 
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Dividing  by  r  and  rearranging  gives  the  equation  of  motion  for  <j)  . 
(44)  <j)  =  —  (  — 2n|)  +  cisina  cos  8  )  -9 2  sintf)  cos(|) 
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APPENDIX  D.  TRANSFORMATION  INTO  HELIOCENTRIC  - 
ECLIPTIC  COORDINATES  (E-FRAME) 


The  matrix  derivations  below  are  used  to  transform  the  spacecraft  state  and  event 
conditions  into  a  common  coordinate  system  -  heliocentric  ecliptic  Cartesian  coordinates, 
called  the  E-frame.  Spacecraft  states  at  the  events  (position  and  velocity)  are  given  in 
spherical  coordinates  while  the  target  planet  manifolds  (elliptical  paths  with  velocity 
related  to  position)  are  most  easily  expressed  in  a  perifocal  frame. 


z 


Figure  57  Spherical  coordinate  system 

SPACECRAFT 

Referring  to  the  Figure  57,  we  transform  the  spherical  position  coordinates  into 
Cartesian  coordinates  as  follows.  It  is  worth  noting  that  some  texts  (e.g.  ref  18)  choose 
the  <j)  coordinate  as  the  complement  angle  from  the  one  depicted  above,  so  it  is  important 
to  be  careful  in  deriving  the  transformation  matrices. 


X 

r  cos(|)  cos  9 

y 

= 

r  cos<|)  sin  9 

|_  z 

E 

1 

-©- 

.5 

c/3 
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The  spacecraft  velocity  state  has  elements  given  as  r,0 ,cf) .  In  the  parlance  of 
DIDO,  these  spacecraft  states  are  the  “primal. states”.  First  we  need  the  velocity  vector  of 
the  spacecraft  in  vr,ve,v1)1  components  where  the  subscripts  indicate  which  sense 
direction  the  velocity  vectors  point.  This  transformation  appears  as 


K 

r 

ve 

= 

?0  cost]) 

3. 

rt |) 

or  in  our  notation 


(45) 


vr 

V,. 

ve 

= 

rcoe  cost]) 

3. 

'Wo 

These  coordinates  may  now  be  transformed  into  the  E-frame  using  a  3-rotation 
and  2-rotation  matrix. 


=  [R3  (-9  )][R2  (4> )] 


We  choose  to  call  the  matrix  [R3  (-9 )]  [  R  ,(<}>) ]  the  B  matrix. 


COS0  COS(|) 

-sin9 

cos0  sin  (() 

B  = 

sin9  cost)) 

COS0 

sin  0  sin  <|> 

—  sin  (f> 

0 

cost]) 

Using  transformations  given  in  boxes,  we  can  transform  the  position  and  velocity 
of  the  spacecraft  in  spherical  coordinates  (primal. states)  into  the  E-frame  coordinate 
frame. 
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PLANETS 


To  define  the  manifolds  of  Earth  and  Mars  we  need  their  positions,  defining  an 
ellipse,  and  equations  for  their  velocities  at  any  point  on  the  ellipse.  Time  is  imposed  as  a 
non-linear  constraint  in  the  event  conditions  without  need  for  transformation,  so  it  will 
not  be  addressed  here.  An  ellipse  may  be  represented  easily  in  either  polar  or  Cartesian 
form,  so  we  choose  the  polar  form  since  most  astrodynamic  textbooks  use  this  form. 
Positions  and  velocities  in  the  perifocal  frame  (Figure  58)  are  related  by  the  following 
equations. 


r  =  ■ 


1  +  <?cos0 . 


and  V  =  m 


2  n 


The  0  _  used  in  these  equations  is  the  angular  positions  of  the  planet  from  their 


perihelion  at  a  particular  event.  The  relation  between  the  planet  true  anomaly  and  the 
spacecraft  angular  position  is 


for  small  inclinations  (ref  4).  For  larger  inclinations  another  transformation  is  required. 

Taking  the  locus  of  planet  positions,  the  ellipse,  and  transforming  in  to  Cartesian 
coordinates  we  write 


xp  =  r cosQp 
y  =rsin0 

J  P  P 

zp=  0 

Transforming  into  the  E- frame  requires  a  3-1-3  rotation  matrix  defined  using  the 
longitude  of  ascending  node  Cl ,  inclination  i ,  and  argument  of  periapsis  to  (ref  13). 

This  matrix  that  transforms  Cartesian  coordinates  from  the  perifocal  frame  to  the 
E- frame  is  defined  as 

£A  =[r3(-co)][^  (-0][R3(-O)] 

cos  Cl  cos  co -cos  (sin  Osin  to  -cos  i  sin  Qcosco  -cos  Osin  co  sinisinco 
=  cos  i  cos  Cl  sin  co  +  sin  Qcosco  cos  i  cos  Qcosco-  sin  Cl  sin  co  -sin  /cos  Cl 
sin  i  sin  co  sin  i  cos  co  cos  i 
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R  denotes  rotation  matrices  about  the  1,  2,  or  3  axis.  See  ref  13,  p.80,  ref  18, 
p.190,  and  ref  4,  p.  151  for  more  details. 

Now  transforming  the  planetary  velocity  vectors  from  the  perifoc  al  coordinates 
(Figure  58),  u  and  v  in  figure  to  the  E- frame.  First,  transform  u  and  v  into  Cartesian 
coordinates  in  the  perifocal  frame. 


V, 

U 

V 

=  [~R,(-e  )1 

V 

y 

[_  3  V  P  >  J 

V 

0 

L  z  J 

perifocal 

Finally,  transforming  the  Cartesian  velocity  vector  in  the  perifocal  plane  into  the 
E- frame,  we  use  the  A  matrix  again. 
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APPENDIX  E.  SOLAR  SAIL  CONTROL  HODOGRAPH  ANALYSIS 


A  hodograph  is  useful  in  determining  convergence  char  acteristics  of  a  solution 
using  a  particular  set  of  controls  in  the  dynamics.  Figure  59  shows  the  hodograph,  the 
mapping  of  the  control,  a,  onto  the  u-v  plane.  Because  the  thrust  magnitude  and 
resulting  acceleration  due  to  the  sail  is  completely  coupled  with  the  sail  angle,  the 
controls  map  as  closed  curve,  not  a  region.  Equations  (38)  and  (40)  plotted  for  all 
feasible  0C  values  produces  the  hodograph.  The  curve  is  not  convex  in  that  a  line  drawn 
between  any  two  points  on  the  curve  does  not  remain  completely  in  the  locus  of  points 
defined  by  the  control  mapping.  There  may  be  a  relation  of  convexity  and  convergence 
speed. 


V 


Figure  59  Hodograph  for  Cone  Angle,  a 

Using  a  new  control  scheme  by  setting  the  control  variable  equal  to  a  maps  as  a 
straight  line,  which  is  inherently  convex  using  the  above  definition.  Solutions  do  in  fact 
converge  faster  for  this  case. 
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